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CHAPTER  I 


Introduction  and  Motivation 

1.1  Introduction 

Hypersonic  flight  vehicles  are  a  current  topic  of  interest  in  both  civilian  and  mili¬ 
tary  research.  NASA  is  currently  designing  a  Crew  Transport  Vehicle  (CTV)  [44,  69] 
and  Crew  Exploration  Vehicle  (CEV)  [32]  to  replace  the  space  shuttle;  reentry  vehi¬ 
cles  are,  by  definition,  hypersonic  vehicles.  Military  requirements  for  reconnaissance 
and  surveillance,  as  well  as  the  mission  of  the  United  States  Air  Force  to  rapidly 
project  power  globally  makes  the  design  of  a  hypersonic  plane  that  can  quickly  tra¬ 
verse  the  globe  very  attractive  [102], 

The  design  of  hypersonic  vehicles  requires  accurate  prediction  of  the  surface  prop¬ 
erties  while  in  flight.  These  quantities  are  typically  the  heat  flux,  pressure  and  shear 
stress,  from  which  the  aerodynamic  forces  and  moments  can  be  calculated.  These 
variables  govern  not  only  the  aerodynamic  performance  of  the  vehicle,  but  also  deter¬ 
mine  the  selection  and  sizing  of  the  thermal  protection  system  (TPS),  which  protects 
the  vehicle  from  the  extreme  temperatures  encountered  at  hypersonic  velocities. 

The  geometry  of  a  vehicle,  and  in  particular,  the  nose  and  the  leading  edges 
of  wings  and  other  aerodynamic  surfaces,  is  a  critical  consideration  in  a  vehicle’s 
design.  Aerodynamic  heating  is  inversely  proportional  to  the  square  root  of  the 
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radius  at  the  stagnation  point;  hence,  historically  most  vehicles  have  had  blunted 
noses  and  leading  edges  to  reduce  the  thermal  loads  to  acceptable  levels. 

Recently,  however,  a  class  of  materials,  designated  Ultra-High  Temperature  Ce¬ 
ramic  (UHTC)  composites,  has  been  developed  that  can  withstand  temperatures  as 
high  as  3500  K  [57,  78].  Materials  such  as  these  allow  the  use  of  much  sharper  leading 
edges.  Sharp  leading  edges  are  important  in  the  design  of  waveriders,  a  class  of  hy¬ 
personic  vehicles  that  depend  on  the  high  pressures  behind  a  shock  wave  to  achieve 
a  high  lift-to-drag  ratio  [2],  These  vehicles  are  designed  theoretically  with  infinitely 
sharp  leading  edges  in  order  that  the  shock  stays  attached.  Manufacturability  and 
thermal  considerations  then  require  a  finite  amount  of  blunting.  Any  blunting  will 
detach  the  shock  allowing  spillage  of  high  pressure  gases  around  the  leading  edge, 
decreasing  aerodynamic  performance  by  as  much  as  20%  [27].  Other  vehicle  designs, 
such  as  the  experimental  X-43A,  also  depend  on  sharp  leading  edges  [97]. 

During  its  trajectory  through  an  atmosphere,  a  hypersonic  vehicle  will  experience 
vastly  different  flow  regimes  because  the  atmosphere’s  density  varies  as  a  function  of 
altitude.  Flight  testing  and  reproduction  of  these  varied  flow  conditions  in  ground- 
based  laboratory  facilities  is  both  expensive  and  technically  challenging.  Hence, 
there  is  an  extremely  important  role  for  computational  models  in  the  development 
of  hypersonic  vehicles. 

1.2  Nonequilibrium  Hypersonic  Gas  Flows 

There  are,  generally  speaking,  three  regimes  in  which  hypersonic  vehicles  travel. 
They  are  classified  as  the  continuum,  continuum-transition  and  free-molecular  re¬ 
gimes.  Typically,  the  different  regimes  are  distinguished  by  the  Knudsen  number, 
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Figure  1.1:  The  Knudsen  number  limits  on  the  mathematical  models  (from  Ref.  [14], 
after  Ref.  [4]). 


defined  as 

(1-1) 

where  the  mean  free  path,  A,  of  a  gas  is  defined  as  the  average  distance  a  particle 
travels  between  successive  collisions;  L  is  a  characteristic  length  and  p  is  the  density. 

Figure  1.1  illustrates  the  different  flow  regimes,  from  the  continuum  (Kn  <  0.01), 
through  the  transitional  to  the  free-molecular  (Kn  >  10)  regime.  The  Boltzmann 
equation  is  valid  for  all  regimes,  whereas  the  Navier-Stokes  equations  are  valid  only 
for  those  regimes  near  the  continuum  limit.  Extended  hydrodynamic  equations,  also 
known  as  higher-order  moment  equations,  can  be  used  further  into  the  transitional 
regime  than  the  Navier-Stokes  equations.  Methods  based  on  these  types  of  equations, 
however,  are  not  as  mature  and  are  subject  to  significant  limitations  that  prevent 
their  use  for  hypersonic  flows,  as  will  be  discussed  in  Chapter  II. 

At  low  altitudes,  the  atmospheric  density  is  relatively  high,  and  flows  around 
hypersonic  vehicles  should  be  simulated  using  traditional  Computational  Fluid  Dy- 
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namics  (CFD)  by  solving  either  the  Euler  or  preferably  the  Navier-Stokes  equations. 
This  is  the  continuum  regime  characterized  by  very  large  Reynolds  numbers  and  very 
low  Knudsen  numbers.  (Although  CFD  refers  to  techniques  used  to  solve  any  set  of 
conservation  equations,  the  term  will  be  used  herein  to  refer  only  to  methods  used 
to  solve  the  Navier-Stokes  equations.) 

At  very  high  altitudes,  at  the  edge  of  the  atmosphere,  the  density  is  low  such  that 
there  are  very  few  collisions  between  the  molecules  and  atoms  in  the  flow  around  the 
vehicle.  This  is  the  rarefied  flow  regime  and  can  be  computed  using  the  direct 
simulation  Monte  Carlo  (DSMC)  method  [4],  which  has  been  shown  to  converge  to 
solutions  of  the  Boltzmann  equation  [95].  Generally  speaking,  CFD  methods  are 
about  an  order  of  magnitude  faster  than  the  DSMC  method  (although  the  DSMC 
method’s  computational  cost  decreases  in  more  rarefied  flows).  However,  the  lack  of 
collisions  makes  the  physics  of  the  Navier-Stokes  equations  invalid  in  rarefied  regimes, 
which  are  characterized  by  a  large  Knudsen  number.  On  a  blunt  body,  a  high-density 
fore-body  flow  can  create  a  rarefied  flow  in  the  wake  of  the  vehicle.  In  principle,  the 
DSMC  method  can  be  applied  to  any  dilute  gas  flow,  but  becomes  prohibitively 
expensive  for  Knudsen  numbers  less  than  0.001.  Thus,  it  is  attractive  to  find  ways 
to  increase  the  validity  of  CFD  methods  beyond  the  continuum  regime. 

One  way  to  improve  CFD  modeling  in  the  transition  regime,  that  is,  for  lower 
density  flows  beyond  the  continuum  regime  and  before  the  free  molecular  regime,  is 
by  replacing  the  typical  no-slip  boundary  conditions  with  slip  velocity  and  temper¬ 
ature  jump  boundary  conditions.  The  addition  of  slip  boundary  conditions  will  not, 
however,  eliminate  all  source  of  errors  when  using  continuum  methods  for  flows  with 
large  amounts  of  nonequilibrium. 

Hybrid  methods  in  which  the  computational  domain  is  split  between  particle 
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(DSMC)  and  continuum  (CFD)  methods  are  another  way  to  decrease  computational 
cost  while  maintaining  accuracy  and  are  an  area  of  active  research  [82], 

The  areas  of  the  flow  where  the  continuum  hypothesis  breaks  down  (or  equiv¬ 
alently,  where  the  flow  is  no  longer  in  local  thermodynamic  equilibrium),  can  be 
quantified  by  the  use  of  a  continuum  breakdown  parameter  [11], 

1.3  Survey  of  Recent  and  Current  Research 

The  DSMC  and  CFD  methods  have  both  been  used  for  many  years  to  model  hy¬ 
personic,  nonequilibrium  gas  flows,  while  advances  in  computing  technology  during 
the  past  several  years,  as  well  as  more  sophisticated  models  have  enabled  the  mod¬ 
eling  of  more  complex  nonequilibrium  phenomena  [8,  42,  67,  70,  91,  104,  106,  108]. 

Some  of  the  earliest  work  comparing  DSMC  and  Navier-Stokes  simulations  is  that 
of  Moss  and  Bird  [60] ,  originally  reported  in  1984,  in  which  DSMC  solutions  of  the 
shuttle  orbiter  nose  during  re-entry  were  compared  with  viscous  shock  layer  (VSL) 
solutions.  Their  work  included  a  5-species  chemistry  model,  and  the  results  showed 
that  there  was  reasonably  good  agreement  in  flow  solutions  at  the  lower  altitudes, 
which  worsened  as  the  flow  became  more  rarefied. 

The  1990’s  saw  an  increase  in  the  use  of  DSMC,  as  well  as  additional  comparisons 
between  continuum  and  particle  solvers.  Only  a  representative  sample  of  this  work 
is  described  here. 

Moss,  et  al.  [59,  62]  compared  DSMC  and  CFD  solutions  of  a  Mach  20  flow 
of  non-reacting  and  reacting  nitrogen  about  a  70-deg  blunted  cone,  for  freestream 
Knudsen  numbers  of  0.001,  0.01  and  0.03.  An  emphasis  was  placed  on  the  wake 
structure  and  afterbody  heating.  A  three  temperature  model  and  slip  boundary  and 
no-slip  boundary  conditions  were  used  for  the  Navier-Stokes  solutions. 
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Olynick,  et  al.  [64]  compared  DSMC  and  CFD  solutions  of  the  flow  about  the 
Fire  II  experimental  re-entry  vehicle.  Freestream  flow  Knudsen  numbers  considered 
were  about  0.0025  and  0.01.  The  solutions  used  a  5-species  model  for  reacting  air. 
Separate  translational,  rotational  and  vibrational  energy  equations,  along  with  slip 
boundary  conditions,  were  employed  in  the  CFD  computations.  Special  attention 
was  paid  to  the  submodels  used  in  both  DSMC  and  CFD  to  ensure  compatibility. 
Emphasis  was  on  the  flow  field  solutions. 

Candler,  et  al.  [16]  simulated  the  flow  of  a  spherical  blunt-body  during  re-entry 
using  CFD  and  DSMC  methods  and  compared  the  resulting  radiative  emissions  from 
the  flow  field.  Similar  chemical  kinetics  and  thermal  relaxation  models  were  used, 
where  possible. 

Research  into  developing  hybrid  DSMC-CFD  methods  also  prompted  additional 
comparisons  between  CFD  and  DSMC  solutions.  Boyd,  et  al.  [11]  computed  so¬ 
lutions  about  a  blunt  sphere  and  for  ID  shockwaves  while  Hash  and  Hassan  [35] 
computed  the  flow  about  a  70-deg  blunted  cone. 

Carlson,  et  al.  [20]  conducted  CFD  and  DSMC  simulations  of  a  hemisphere  in 
air  for  Mach  10  and  15.  The  flow  was  very  near  continuum,  with  Knudsen  numbers 
of  0.02  and  below.  Different  models  of  air  chemistry  (perfect  gas,  equilibrium  air,  5- 
species)  were  used,  as  well  as  some  vibrational  nonequilibrium  models.  The  emphasis 
was  on  the  effect  of  thermochemical  nonequilibrium  on  the  field  of  view  of  a  sensor. 

More  recently,  work  on  comparing  DSMC  and  CFD  solutions  has  been  con¬ 
centrated  on  several  validation  cases  for  computational  codes  [107].  In  particular, 
the  NATO  Research  Technology  Organization  (RTO)  Advanced  Vehicle  Technology 
Panel  Working  Group  10  (WG  10)  coordinated  several  experiments  to  highlight  six 
topical  areas  for  CFD  validation  [96].  These  areas  included  shock-shock  interactions 
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and  laminar  hypersonic  viscous-inviscid  interactions  [45].  Computational  simula¬ 
tions  for  two  experiments,  Mach  12  and  Mach  16  flows  of  nitrogen  over  a  25-55-deg 
double  cone  and  a  hollow  cylinder  flare,  were  solicited  for  a  blind  comparison  for  the 
2001  AIAA  Aerospace  Sciences  Meeting  and  Exhibit.  Several  DSMC  [13,  58]  and 
CFD  [18,  28,  43]  solutions  were  submitted. 

Later,  others  compared  CFD  and  DSMC  solutions  to  these  and  similar  experi¬ 
mental  results  [19,  55,  61,  72,  73,  100,  99].  Inger  and  Moss  [39]  also  compared  DSMC 
with  theoretically  derived  expressions  from  the  Navier-Stokes  equations  for  the  sepa¬ 
ration  and  reattachment  streamline  angles  for  the  shock-boundary  layer  interaction. 

In  each  of  these  validation  cases,  while  the  surface  properties  between  CFD  and 
DSMC  were  computed  and  compared  with  the  experimental  results,  particular  em¬ 
phasis  was  on  the  size  of  the  recirculation  zone  near  the  shock-boundary  layer  inter¬ 
action. 

Most  recently,  Boyd,  et  al.  [12]  and  Ozawa,  et  al.  [65]  compared  particle  and 
continuum  solutions  with  flight  data  for  the  Stardust  atmospheric  re-entry  for  near¬ 
continuum  conditions  (with  a  freestream  Knudsen  number  of  about  0.005).  This 
data  set  is  of  particular  interest  considering  the  high  velocities  (about  12.6  km/s)  at¬ 
tained  during  re-entry.  This  study  focused  on  dissociation  and  ionization.  Enormous 
differences  were  seen  in  basic  flow  property  predictions  between  the  two  methods. 

Jain  and  Hayes  [40]  developed  an  analytical  method  for  engineering  estimates  of 
pressure,  shear  stress  and  heat  transfer  rates  on  vehicles  of  arbitrary  shape  for  the 
hypersonic  continuum  through  the  transitional  and  free-molecular  flow  regimes.  The 
method  is  applicable  to  sharp-  and  blunt-nosed  bodies.  Solutions  are  compared  with 
DSMC  and  CFD  solutions  with  reasonable  accuracy. 

There  is,  then,  an  abundant  amount  of  research  that  has  been  performed,  and 


is  still  being  done,  to  determine  how  accurate  CFD  solutions  are  compared  with 
DSMC  and  experimental  data.  However,  these  studies  are  limited  in  several  ways; 
for  example,  few  studied  cases  for  several  different  flow  regimes  or  body  geometry.  In 
addition,  as  most  are  compared  to  experiment,  they  include  complex  thermochemical 
nonequilibrium  effects.  The  addition  of  these  complex  models,  while  important,  can 
mask  fundamental  differences  that  must  be  understood.  Furthermore,  while  the 
surface  properties  are  computed  in  some  studies,  the  emphasis  is  usually  placed 
on  other  flow  properties,  such  as  ionization  species  concentration  and  shock  wave 
structure. 

A  quantitative  link  between  a  given  level  of  continuum-breakdown  and  the  ac¬ 
curacy  of  predicted  surface  quantities  using  CFD  has  not  been  presented  in  prior 
studies.  Thus,  there  is  a  need  for  a  more  systematic,  fundamental  study  to  deter¬ 
mine  the  effects  of  nonequilibrium  on  the  surface  properties  of  hypersonic  vehicles. 
The  goal  of  the  present  study  is  therefore  to  investigate  this  issue.  Specifically, 
how  are  the  critical  hypersonic  vehicle  design  surface  properties  of  pressure,  shear 
stress  and  heat  transfer  rate  affected  by  failure  of  the  continuum  approach  in  certain 
regions  of  the  flow  field?  For  example,  in  hypersonic  flow,  the  first  place  where  con¬ 
tinuum  breakdown  is  observed  is  within  the  shock  wave  itself.  It  is  well  known  that 
traditional,  Navier-Stokes-based  CFD  cannot  accurately  predict  hypersonic  shock 
structure  [17,  25].  It  is  not  clear,  however,  whether  local  breakdown  within  the 
shock  has  a  tangible  impact  on  the  rest  of  the  flow  field  and  the  resulting  surface 


properties. 
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1.4  Scope  of  Current  Work 

The  research  presented  in  this  dissertation  has  several  goals  meant  to  address 
the  limitations  in  the  previous  and  current  research  efforts  mentioned  above.  These 
goals  are,  specifically: 

1.  Start  with  the  fundamentals.  The  present  study,  as  a  purely  numerical  study, 
will  focus  primarily  on  the  fundamentals  of  nonequilibrium  behavior  and  grad¬ 
ually  increase  the  complexity,  starting  with  a  monatomic  gas,  argon,  and  pro¬ 
gressing  to  a  diatomic  gas,  nitrogen.  The  effects  of  each  type  of  nonequilibrium 
on  the  surface  properties  will  be  quantified  as  the  complexity  increases. 

2.  Study  many  flow  regimes,  about  blunt  and  sharp  bodies.  The  current  work 
will  consider  flow  regimes  from  the  continuum  and  into  the  transitional  regime 
to  quantify  the  effects  of  the  degree  of  rarefaction;  considering  two  different 
flow  velocities  to  quantify  the  effects  of  larger  Mach  number;  and  considering 
two  types  of  geometry,  a  cylinder  and  a  wedge,  to  quantify  differences  clue  to 
blunt-body  phenomena  versus  sharp  leading-edge  phenomena. 

3.  Evaluate  the  effectiveness  of  several  types  of  CFD  slip  boundary  conditions 
and  compare  the  CFD  slip  values  with  the  DSMC  slip  values.  This  research 
will  evaluate  the  effectiveness  of  several  CFD  slip  boundary  conditions,  includ¬ 
ing  one  only  recently  proposed  [49],  in  predicting  the  surface  properties  of  a 
hypersonic  vehicle.  The  actual  slip  quantities  predicted  by  these  boundary  con¬ 
ditions  will  also  be  compared  with  those  extracted  from  the  DSMC  simulations 


for  each  flow  condition. 
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4.  Lay  the  foundation  for  further  studies  essential  to  the  design  of  hybrid  methods. 
Hybrid  methods  face  two  basic  problems;  determining  the  boundaries  between 
the  CFD  and  DSMC  domains  and  passing  information  from  one  domain  to 
the  other.  This  research  contributes  to  both  of  these  areas.  The  chosen  value 
for  the  continuum  breakdown  parameter’s  effectiveness  in  predicting  differences 
will  be  shown  by  comparing  the  breakdown  parameter  value  with  the  other  flow 
properties.  An  effective  hybrid  design  also  requires  that  the  different  submodels 
used  in  both  computational  methods  be  equivalent  as  much  as  possible;  thus 
information  passed  between  both  domains  is  also  equivalent. 

5.  Show  conclusively  that  flow  property  differences  near  the  wall  are  concentrated 
in  the  Knudsen  layer.  Unique  to  this  dissertation  are  the  results  that  the 
differences  between  CFD  and  DSMC  near  the  wall  are  concentrated  mainly  in 
the  Knudsen  layer,  defined  here  as  the  region  of  flow  10  mean  free  paths  or  less 
from  the  wall  surface. 

An  outline  of  this  dissertation  is  as  follows: 

Chapter  2  presents  a  brief  description  of  kinetic  theory  and  the  concepts  of  equi¬ 
librium  and  nonequilibrium  gas  dynamics  and  the  equations  governing  gas  flows. 
The  chapter  concludes  with  a  brief  description  of  the  DSMC  code  MONACO  and 
the  CFD  code  LeMANS,  which  are  used  for  the  computational  analyses  in  the  re¬ 
maining  chapters. 

Chapter  3  discusses  the  different  submodels  present  in  DSMC  and  CFD  simu¬ 
lations.  Such  physical  models  include  transport  properties  (such  as  viscosity),  wall 
boundary  conditions  and  vibrational  relaxation.  This  chapter  will  discuss  the  rele¬ 
vant  physical  models  and  the  manner  in  which  they  are  treated  in  each  simulation 
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method  such  that  they  are  equivalent,  as  much  as  is  possible. 

Chapter  4  presents  solutions  obtained  using  both  computational  methods  for  a 
hypersonic  flow  about  a  cylinder.  First,  the  case  of  a  hypersonic  flow  of  argon,  a 
monatomic  gas,  is  considered.  The  monatomic  nature  of  argon  eliminates  the  possi¬ 
bility  of  thermal  nonequilibrium  due  to  the  nonexistence  of  internal  energy  modes. 
Then  the  case  of  a  hypersonic  flow  of  nitrogen,  using  the  same  physical  geometry,  is 
considered.  The  use  of  nitrogen,  a  diatomic  gas,  allows  the  investigation  of  the  effects 
of  thermal  nonequilibrium  in  addition  to  the  translational  nonequilibrium  present  in 
the  argon  flow. 

Chapter  5  considers  the  flow  about  a  wedge  with  a  sharp  leading-edge.  Again, 
solutions  for  flows  of  argon  and  nitrogen  at  Mach  10  and  Mach  25  are  computed  and 
compared.  The  distinct  physical  phenomena  associated  with  a  sharp  leading-edge 
flow  are  discussed. 

Chapter  6  presents  two-dimensional  CFD  solutions  that  are  compared  with  ex¬ 
perimental  measurements  of  a  hypersonic  flow  of  nitrogen  over  a  flat  plate.  Several 
different  values  for  the  accommodation  coefficient  are  evaluated.  In  addition,  the 
CFD  solutions  are  also  indirectly  compared  to  DSMC  solutions  of  the  same  flow. 
Thus,  the  relative  accuracy  of  CFD  and  DSMC  can  be  evaluated  for  a  realistic  flow. 

The  dissertation  concludes  with  Chapter  7  in  which  some  conclusions  are  drawn 


and  future  work  is  proposed. 


CHAPTER  II 


Simulation  of  Hypersonic  Gas  Flows:  Background 

and  Theory 


2.1  Introduction 

The  computational  simulation  of  nonequilibrium  hypersonic  gas  flows  requires  a 
basic  understanding  of  the  kinetic  theory  of  gases,  as  well  as  the  different  methods 
used  to  model  the  varied  phenomena  present.  This  chapter  begins  by  presenting  a 
brief  description  of  kinetic  theory  and  the  concepts  of  equilibrium  and  nonequilibrium 
gas  dynamics.  It  then  describes  the  equations  governing  gas  flows,  including  the 
Boltzmann  equation  and  the  Navier-Stokes  equations.  The  Boltzmann  equation 
can  be  emulated  using  the  direct  simulation  Monte  Carlo  (DSMC)  method,  while 
the  Navier-Stokes  equations  can  be  solved  numerically  using  Computational  Fluid 
Dynamics  (CFD)  techniques.  The  chapter  concludes  with  a  brief  description  of 
the  DSMC  code  MONACO  and  the  CFD  code  LeMANS,  which  are  used  for  the 
computational  analyses  in  the  remaining  chapters. 

2.2  Some  Basics  of  Kinetic  Theory 

The  discussion  of  the  methods  to  follow  requires  some  degree  of  knowledge  in 
the  kinetic  theory  of  gases.  This  section  is  only  meant  to  provide  a  short  overview. 
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Monatomic  gases  (such  as  argon)  and  diatomic  gases  (such  as  nitrogen)  are  consid¬ 
ered.  For  more  details,  the  reader  is  directed  to  many  texts  on  the  subject,  such  as 
References  [4,  31,  93]. 

Kinetic  theory  considers  a  gas  flow  on  the  molecular  level.  The  individual  gas 
molecules,  or  particles,  are  considered  to  be  constantly  moving  about,  colliding  with 
other  particles  and  any  surfaces  present.  The  properties  of  the  flow  depend  only  on 
the  mass,  size,  position,  velocity  and  internal  energy  of  the  particles.  In  this  chapter, 
the  mass  of  a  particle  is  m;  the  size  is  defined  as  an  effective  particle  diameter,  d\ 
the  position  is  given  as  a  vector  position  from  the  origin,  and  the  velocity  of 
an  individual  gas  particle  is  denoted  as  a  vector,  Cj.  The  macroscopic  thermody¬ 
namic  properties,  such  as  temperature,  density  and  pressure,  are  derived  by  taking 
moments,  or  averages,  of  the  individual  particle  properties. 

The  individual  molecular  velocity  can  be  split  into  its  random  and  average  com¬ 
ponents  as  c'  =  | Ci  —  {a)\  where  (c*)  is  the  average  velocity  of  the  particles  in  the 
volume  under  consideration.  The  random  velocity,  c',  is  also  known  as  the  thermal 
velocity  and  the  average  velocity  is  known  as  the  bulk  velocity. 

Each  particle  may  have  several  energy  modes.  The  translational  energy  is  de¬ 
scribed  by  the  random  motion  of  the  particles.  Diatomic  particles  also  possess  inter¬ 
nal  energy  due  to  rotation  of  the  atoms  around  an  axis,  as  well  as  vibration  of  the 
atoms  along  the  internuclear  axis. 

The  thermodynamic,  or  translational,  temperature  can  be  defined  as  a  measure 
of  the  kinetic  energy  due  to  the  random  motion  of  the  gas  particles  and  is  defined  as 

etra  =  \rn  ((d)  +  (c'l)  +  (d))  =  ^k Ttra  (2.1) 


where  etra  is  the  average  translational  energy  per  particle  and  k  is  the  Boltzmann  con- 
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stant.  Similarly,  rotational  and  vibrational  temperatures  can  be  defined  as  measures 
of  the  internal  rotational  and  vibrational  energy  of  a  diatomic  gas 

&rot  =  k  Trot,  (2-2) 

_  k©„ 

evib  —  777  777  \  Tj  'A-'H 

exp  y\3v/Tvn)j  1 

where  erot  and  em^  are  the  average  rotational  and  vibrational  energies  per  particle, 

and  0,,  is  the  characteristic  vibrational  temperature  (the  temperature  at  which  the 
vibrational  mode  is  significantly  activated,  approximately  3300  K  for  nitrogen).  Here 
the  assumption  has  been  made  that  the  rotational  mode  is  fully  activated  at  the 
temperatures  of  interest  (the  characteristic  temperature  of  rotation  for  nitrogen  is 
about  3  K),  and  that  the  vibrational  energy  can  be  modeled  as  a  harmonic  oscillator. 

The  rotational  and  vibrational  energy  modes  are  activated  through  the  process 
of  intermolecular  collisions.  As  the  molecules  collide,  energy  is  transferred  from  the 
translational  mode  to  the  rotational  and  vibrational  modes,  and  vice-versa.  The 
number  of  collisions  required  to  activate  the  internal  energy  modes  is  dependent  on 
the  temperature  of  the  flow.  As  the  temperature  increases,  the  collisions  tend  to  be 
more  energetic,  and  hence,  the  rotational  and  vibrational  modes  are  activated  with 
fewer  collisions. 

Figure  2.1  illustrates  how  the  rotational  and  vibrational  collision  probabilities 
for  nitrogen  vary  with  temperature.  Here,  the  rotational  collision  probabilities  are 
obtained  from  Lordi  and  Mates  experimental  values  [4]  and  the  vibrational  colli¬ 
sion  probabilities  plotted  are  obtained  from  the  Landau-Teller  relaxation  model  as 
explained  in  Chapter  III.  The  collision  probability  is  the  inverse  of  the  collision 
number,  or  the  number  of  collisions  required  on  average  to  activate  the  internal  en¬ 
ergy  modes.  At  lower  temperatures — on  the  order  of  100  K — the  rotational  collision 
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Rotation 


Figure  2.1:  Probability  that  a  collision  will  results  in  transfer  of  rotational  and  vi¬ 
brational  energy  (for  nitrogen). 


probability  is  about  1/3.  That  is,  about  3  collisions  are  required,  on  average,  to 
activate  the  rotational  modes.  As  the  temperature  increases,  the  rotational  colli¬ 
sion  probability  decreases  and  eventually  remains  constant  at  a  value  of  about  0.045 
(equivalent  to  a  collision  number  of  about  21). 

Alternatively,  the  probability  of  vibrational  energy  exchange  remains  very  small 
(below  10-3,  or  a  collision  number  of  about  1000)  until  the  temperature  reaches 
about  5000  K.  As  the  temperature  increases,  the  vibrational  collision  probability 
increases  dramatically  and  levels  off  at  about  0.025  at  50,000  K.  Thus,  at  such  ele¬ 
vated  temperatures,  a  larger  number  of  collisions  would  result  in  vibrational  energy 
exchange  than  at  lower  temperatures.  It  is  important  to  note  that  the  number  of 
collisions  required  for  vibrational  activation  decreases  dramatically  as  the  tempera¬ 
ture  increases,  but  the  rotational  modes  still  require  far  fewer  collisions  for  activation 
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even  at  the  higher  temperatures. 

An  exact  representation  of  a  collision,  or  interaction  of  two  or  more  particles, 
requires  a  detailed  knowledge  of  the  shape  and  orientation  of  the  individual  particles. 
For  any  realistic  gas  flow,  this  is  impossible.  Different  models  for  the  shape  of  the 
intermolecular  force  simplify  the  analysis.  The  hard-sphere  model  describes  each  gas 
particle  as  an  elastic  sphere  with  a  specific  size,  defined  as  the  diameter,  d.  There 
is  no  intermolecular  force  until  the  two  molecules  come  into  contact,  at  which  point 
the  repulsive  force  is  infinite. 

A  significant  weakness  of  the  hard-sphere  model  is  the  fixed  size,  d.  The  total 
collision  cross-section,  o>r,  for  the  hard-sphere  model  is  given  by  cr?  —  n d 2 .  Expe¬ 
rience  has  shown  that  the  total  collision  cross-section  is  dependent  on  the  relative 
speed  between  the  molecules  involved  in  the  collision  and  it  is  important  to  reproduce 
this  behavior  to  successfully  model  the  temperature  dependence  of  viscosity  [4],  The 
average  relative  speed  is  dependent  on  temperature;  hence,  there  is  a  temperature 
dependence  to  the  average  total  collision  cross-section,  and  the  particle  diameter. 

Other  models  have  been  proposed  that  model  the  temperature  dependence  of  the 
collision  cross-section  (and  the  viscosity)  in  a  more  realistic  manner.  Among  these  is 
the  variable  hard  sphere  (VHS)  model  [4],  A  VHS  particle  has  a  diameter  that  is  a 
function  of  the  relative  velocity  of  the  collision  partners.  In  many  cases  the  function 
is  an  inverse  power  law,  with  the  temperature  dependence  explicitly  chosen  to  match 
experimental  viscosity  data. 

In  any  realistic  flow  it  is  impossible  to  follow  each  individual  molecule  as  it  collides 
with  other  molecules  and  surfaces  and  describe  its  particular  properties  as  a  function 
of  time.  The  use  of  a  velocity  distribution  function  (VDF)  allows  a  probabilistic 
description  of  a  particle’s  velocity  and  position.  A  VDF,  denoted  by  /  =  /(q,  Xi,  t ),  is 
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simply  a  probability  density  function  in  velocity-space;  that  is,  it  gives  the  probability 
that  any  particular  particle  will  have  a  velocity  that  falls  within  the  range  c*  +  Ac 
at  a  particular  location  and  time. 

Once  the  VDF  for  a  particular  flow  is  known,  the  macroscopic  properties  can  be 
obtained  by  taking  moments  of  the  VDF.  The  moment  of  a  quantity,  Q(ci),  is  defined 
as 

/OO 

Q(ci)f(ci)dci.  (2.4) 

-OO 

When  Q  =  c"  (where  n  is  a  power),  then  the  average  (c”)  is  known  as  the  nth-nioment 
of  the  VDF. 

2.3  Equilibrium  and  Nonequilibrium 


Choosing  a  simulation  method  appropriate  for  a  particular  gas  flow  depends  on 
whether  or  not  there  are  significant  equilibrium  effects  present.  A  precise  definition 
of  thermodynamic  equilibrium  will  not  be  given  here,  but  rather  a  few  qualitative 
descriptions  will  be  given  to  help  in  understanding  the  difference  between  equilibrium 
and  nonequilibrium  gas  flows. 

A  gas  in  equilibrium  can  be  thought  of  as  one  whose  molecular  properties  are 
unchanging  in  time  and  space.  This  suggests  that  there  are  no  gradients  in  molec¬ 
ular  or  macroscopic  properties  (velocity,  temperature,  mass  density,  etc).  A  gas  in 
equilibrium  will  have  a  velocity  distribution  as  given  by  Maxwell, 


(c'i  +  c; 


1  +  c 


(2.5) 


A  gas  flow  that  is  in  complete  thermodynamic  equilibrium  would  have  no  inter¬ 
esting  features  present,  and  would,  in  fact,  be  completely  at  rest.  The  driving  force 
behind  flow  features  of  interest  are  inherently  nonequilibrium.  However,  if  changes 
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in  the  gas  that  are  clue  to  nonequilibrium  effects  occur  significantly  rapidly  such  that 
the  gas  can  be  thought  of  as  adjusting  instantaneously  to  those  changes,  the  flow 
is  said  to  be  in  local  thermodynamic  equilibrium  (LTE).  Thus,  although  the  gas  is 
not  at  rest,  the  departure  from  the  Maxwellian  distribution  2.5  is  everywhere  small. 
The  assumption  of  LTE  implies  that  the  effects  of  viscosity  and  thermal  conductivity 
are  negligible;  that  is,  there  is  no  transport  of  momentum  or  thermal  energy  due  to 
velocity  and  temperature  gradients. 

In  contrast,  a  gas  in  nonequilibrium  will  have  gradients  that  allow  for  the  trans¬ 
port  within  the  gas  of  mass,  momentum  and/or  energy.  Mass  flow  is  driven  by  species 
concentration  gradients.  A  viscous  fluid  with  velocity  gradients  will  cause  a  trans¬ 
fer  of  momentum.  Similarly,  a  gas  with  a  temperature  gradient  will  transport  heat 
energy.  Thermal  nonequilibrium  concerns  internal  energy  modes  of  a  diatomic  gas, 
specifically  the  rotational  and  vibrational  energy  modes.  Chemical  nonequilibrium, 
present  in  a  reacting  flow,  is  beyond  the  scope  of  the  current  research  and  will  not 
be  discussed  further  here.  Mass  transport  due  to  species  concentration  gradients  is 
also  neglected  as  only  simple  gases,  comprising  one  species,  are  considered. 

The  transfer  of  momentum  and  energy  is  due  to  translational  nonequilibrium 
and  gives  rise  to  the  effects  of  viscosity  and  thermal  conductivity.  Rotational  none¬ 
quilibrium  is  manifested  in  two  different  areas:  the  first  concerns  the  so-called  bulk 
viscosity;  and  the  second  a  thermal  nonequilibrium  where  the  rotational  temper¬ 
ature,  as  defined  in  Eq.  2.2,  is  not  the  same  as  the  translational  temperature,  as 
defined  in  Eq.  2.1.  Similarly,  vibrational  nonequilibrium  is  present  in  gas  flows  where 
the  vibrational  temperature,  as  defined  in  Eq.  2.3,  is  not  equal  to  the  translational 
temperature. 

When  nonequilibrium  effects  are  present,  the  gas  is  driven  towards  equilibrium 


19 


through  intermolecular  collisions.  The  number  of  collisions  required  for  gas  molecules 
to  equilibrate  is  dependent  on  the  type  of  nonequilibrium  present.  Translational 
energy  equilibrates  with  only  a  few  collisions.  Rotational  energy  requires  on  the  order 
of  ten  collisions,  while  vibrational  energy  typically  requires  thousands  of  collisions  to 
equilibrate. 

The  relaxation  time  is  the  time  required  for  the  gas  to  come  to  equilibrium. 
The  distance  traveled  during  that  time  by  the  bulk  flow  is  the  relaxation  distance. 
Generally,  the  relaxation  time  is  on  the  order  of  the  mean  collision  time,  and  the 
relaxation  distance  is  on  the  order  of  the  mean  free  path  [4]. 

The  residence  time  of  a  gas  particle  can  be  defined  as  the  amount  of  time  taken 
for  the  gas  particle  to  traverse  a  given  flow  feature,  such  as  a  velocity  or  temperature 
gradient.  If  the  residence  time  is  much  longer  than  the  relaxation  time — that  is,  if 
the  molecules  undergo  sufficient  collisions  to  equilibrate  to  the  local  thermodynamic 
properties — then  the  flow  is  in  equilibrium.  However,  if  the  residence  time  is  short 
compared  to  the  relaxation  time,  the  gas  particle  will  not  reach  equilibrium  with  the 
local  thermodynamic  conditions.  Therefore,  nonequilibrium  effects  can  be  expected 
in  flow  conditions  with  low  residence  times,  or  in  flow  conditions  where  relaxation 
times  are  large.  These  conditions  are  seen  in  areas  of  large  gradients  (such  as  in  a 
shock  wave  or  boundary  layer)  and  in  rarefied  conditions  (where  the  mean  free  path, 
and  mean  collision  time,  is  large). 

In  a  hypersonic  flow,  then,  there  are  three  main  causes  of  significant  nonequilib¬ 
rium 


•  High  velocities  result  in  shorter  residence  times  and  larger  gradients. 


•  High  temperatures  activate  the  vibrational  energy  modes,  which  are  slower  to 
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equilibrate  than  other  energy  modes. 

•  Many  hypersonic  vehicles  fly  in  the  upper,  rarefied  atmosphere.  As  the  density 
decreases,  the  mean  collision  time  increases  and,  thus,  the  relaxation  time 
increases. 

2.4  The  Governing  Equations  of  Gas  Flows 

This  section  describes  the  governing  equations  of  gas  flow,  starting  with  the  Boltz¬ 
mann  equation.  The  manner  in  which  the  Navier-Stokes  equations  can  be  derived 
from  the  Boltzmann  equation  is  then  briefly  described.  The  resulting  discussion 
highlights  the  strengths  and  weaknesses  of  each  approach  when  used  to  model  non¬ 
equilibrium  gas  flows. 

2.4.1  The  Boltzmann  Equation 

The  Boltzmann  equation,  Eq.  2.6,  describes  the  evolution  in  phase  space  (a 
combination  of  velocity  and  physical  space)  of  the  velocity  distribution  function  of 
a  particular  gas  flow  [93].  There  are  two  convective  terms  present;  one  models  the 
convection  in  physical  space  due  to  the  velocity,  cy  and  the  other  the  convection  in 
velocity  space  due  to  accelerations  caused  by  a  force,  Fy  The  source  term  models 
the  increase  and  decrease  of  particles  of  class  c*  due  to  collisions. 

^K(c0]  +  ci^-K(c*)]  +  7^[*>/(ci)]  =  |^[n/(Q)]|  ^  (2.6) 

The  form  of  the  collision  term  depends  on  the  particular  molecular  model  con¬ 
sidered.  For  a  simple  binary  collision  model,  it  can  be  written  as 

|^K(q)]|  =  I  n2[/(c')/(C')  ~  f(ci)f(Ci)\gadndCi,  (2.7) 
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where  a  is  the  differential  collision  cross-section,  g  is  the  relative  velocity  of  the 
two  colliding  particles  (g  —  c*  —  Q  | ) ,  and  dfl  is  the  differential  solid  angle  associated 
with  the  collision.  Two  types  of  collisions  are  considered.  The  first  involves  collisions 
between  particles  of  class  c%  with  particles  of  class  Q.  These  collisions  deplete  the 
number  of  particles  in  class  q.  The  second  type  of  collisions  is  the  inverse  of  the  first; 
that  is,  collisions  between  particles  of  class  c'  and  class  ('■  These  collisions  replenish 
the  number  of  particles  in  class  ct.  The  total  effect  of  these  collisions  on  the  VDF  is 
found  by  integrating  over  all  solid  angles,  and  all  collision  pair  velocities,  Q. 

The  Boltzmann  equation  is  valid  for  all  regimes  of  a  gas  flow,  from  the  continuum 
to  the  rarefied  regime,  although  it  has  been  derived  above  to  consider  only  binary 
collisions  which  would  limit  its  validity  to  dilute  gases.  The  main  challenge  in  us¬ 
ing  the  Boltzmann  equation  for  modeling  gas  flows  is  the  collision  integral.  Even 
assuming  binary  collisions  only,  the  term  is  impossible  (due  to  its  nonlinear  integral 
nature)  to  solve  analytically  and  difficult  to  model  numerically. 

The  Maxwellian  VDF,  given  in  Eq.  2.5  is  a  solution  to  the  Boltzmann  equation 
when  the  collision  integral  term  is  zero,  and  the  flow  is  considered  to  be  in  LTE 
everywhere. 

The  Moment  Equations 

Moments  of  the  velocity  distribution  function  were  defined  in  Eq.  2.4.  Similarly, 
moments  can  be  taken  of  the  Boltzmann  equation  to  give  the  moment  equations,  or 
equations  of  transfer, 

|(„(Q»  +  A(„ {CjQ))  _  =  A[Q],  (2.8) 

where  A[Q]  is  the  moment  of  the  collision  integral  term. 

It  can  be  shown  [93]  that  when  the  moment,  Q(ci),  is  taken  to  be  the  mass, 
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momentum  or  energy  per  particle  (m,  mCi  or  mcf/ 2),  the  change  in  Q  for  the  collision 
partners  must  remain  zero,  and,  thus,  A[Q]  =  0.  Further  simplifying  the  resulting 
set  of  equations  gives  the  Euler  equations.  The  Euler  equations  are  then  equivalent 
to  the  Maxwellian  VDF  when  taking  moments  of  the  Boltzmann  equation.  The  Euler 
equations  are  appropriate  for  modeling  gas  flows  under  the  assumption  of  LTE. 

The  Chapman-Enskog  Expansion 

As  mentioned,  the  Boltzmann  equation  has  an  equilibrium  solution  of  /  =  /o, 
where  /0  is  the  equilibrium,  or  Maxwellian,  distribution  function  given  in  Eq.  2.5. 
Power-series  solutions  can  be  constructed  for  the  Boltzmann  equation.  One  well- 
known  example  is  the  Chapman-Enskog  solution. 

The  Chapman-Enskog  solution  is  obtained  first  by  nondimensionalizing  the  Boltz¬ 
mann  equation  in  terms  of  a  parameter  £.  It  can  be  shown  [93]  that  the  parameter 
£  is  proportional  to  Kn  =  X/L.  Thus,  for  gas  flows  where  Kn  <<  1,  this  parame¬ 
ter  will  be  small.  As  £  approaches  zero,  /  approaches  /0;  this  equation  describes  a 
small  departure  from  equilibrium  (a  perturbation  model).  The  Euler  equations,  the 
Navier-Stokes  equations  and  the  Burnett  equations  result  from  the  Chapman-Enskog 
expansion  of  the  distribution  function  for  small  departures  from  /0. 

As  a  power  series,  the  VDF  can  be  written  as 

/  =  /o(l  +  £0i  +  £202  +  •••) 

where  /  is  the  non-dimensional  VDF. 

The  series  is  then  usually  truncated  after  one,  two  or  three  terms  and  substituted 
back  into  the  Boltzmann  equation,  of  which  moments  are  taken.  The  resulting 
moment  equations  are  the  Euler  equations  (if  only  one  term  is  kept),  the  Navier- 
Stokes  equations  (if  two  terms  are  kept)  and  the  Burnett  equations  (if  three  terms 
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are  kept). 

2.4.2  The  Navier-Stokes  Equations 


The  Navier-Stokes  equations  (defined  here  as  including  the  mass  and  energy  con¬ 
servation  equations  in  addition  to  the  momentum  conservation  equations)  are  typi¬ 
cally  used  to  describe  gas  flows  in  the  continuum  regime.  As  they  can  be  derived  from 
the  Chapman-Enskog  expansion  of  the  Boltzmann  equation  (by  keeping  first-order 
terms),  they  are  valid  only  for  flows  with  small  perturbations  from  equilibrium. 

The  Navier-Stokes  equations  for  a  simple  gas  and  neglecting  body  forces  can  be 
written  as  [93]: 
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where  E  =  e  +  \u.iUi  is  the  total  energy  per  mass  (e  is  the  internal  energy  per  mass, 
which  includes  the  translational,  rotational  and  vibrational  energy),  and  TtJ  and  q 
are  the  shear  stress  tensor  and  heat  flux  vector,  respectively.  (Note  that  iq  =  (q).) 
In  addition  to  Ttl  and  q,  which  will  be  discussed  shortly,  another  equation  is  required 
to  close  the  set.  Typically,  an  equation  of  state,  such  as  the  perfect  gas  law,  is  used. 

The  shear  stress  tensor  and  heat  flux  vectors  arise  due  to  translational  nonequil¬ 
ibrium,  and  can  be  derived  as 
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(2.13) 


where  p  is  the  coefficient  of  viscosity,  ps  is  the  bulk  viscosity  and  k  is  the  coefficient 
of  thermal  conductivity.  For  a  diatomic  gas,  similar  expressions  for  the  heat  flux  clue 
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to  rotational  and  vibrational  energy  are  needed.  Note  that  the  Enler  equations  can 
be  recovered  from  the  Navier-Stokes  equations  above  by  setting  TtJ  and  qi  equal  to 
zero. 

Although  the  effects  of  the  viscosity  and  thermal  conductivity  are  due  to  trans¬ 
lational  nonequilibrium,  the  effect  of  bulk  viscosity  is  due  to  rotational  nonequilib¬ 
rium  [93].  For  monatomic  gases,  then,  it  is  equal  to  zero.  Most  conventional  fluid 
dynamic  analyses  assume  it  is  zero  even  for  diatomic  gases  (Stokes’  hypothesis). 

The  values  of  the  coefficients  of  viscosity  and  thermal  conductivity  for  a  particular 
gas  can  also  be  derived  from  kinetic  theory.  Collision  integral  calculations  are  used 
to  accurately  determine  their  values  across  a  wide  range  of  temperatures  [105]. 

The  Navier-Stokes  equations  are  only  valid  for  the  continuum  regime  (Kn  <  0.01) 
with  the  no-slip  boundary  condition.  Their  validity  can  be  extended  to  Kn  <  0.1 
by  using  slip  boundary  conditions,  but  for  higher  Knudsen  numbers,  they  fail  to 
accurately  predict  the  flow. 

Higher  Order  Moments  and  Extended  Hydrodynamics 

Additional  information  can  be  derived  from  the  Boltzmann  equation  by  retaining 
higher  order  terms  of  the  Chapman-Enskog  expansion.  The  Burnett  equations  result 
from  retaining  the  first  three  terms  of  the  expansion,  and  the  super-Burnett  equa¬ 
tions  result  from  retaining  the  first  four  terms.  While  the  Burnett  equations  can  give 
a  more  accurate  description  of  the  flow  in  nonequilibrium  flows  (such  as  the  interior 
of  shock  waves  [25]),  there  remain  several  significant  hurdles  to  their  practical  im¬ 
plementation.  Some  of  these  include  numerical  stability  and  a  failure  to  satisfy  the 
second-law  of  thermodynamics  [22],  Some  researchers  also  contend  that  the  Burnett 
equations  cannot  be  used  where  the  Navier-Stokes  equations  have  already  failed,  as 
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they  are  also  valid  only  for  Knudsen  numbers  less  than  unity  [33,  34], 

Other  attempts  at  deriving  a  general  hydrodynamic  approach  for  nonequilibrium 
gas  flows  include  higher  moment  methods,  of  which  Grad’s  method  is  one  example. 
In  this  approach,  higher  order  moments  are  taken  of  the  Boltzmann  equation  and 
related  to  the  lower  order  moments.  A  system  of  20  equations  can  be  obtained, 
which  are  then  simplified  to  a  set  of  13  moment  equations  [31].  Sets  of  higher  order 
moment  equations  can  be  obtained  from  the  Boltzmann  equation  by  considering 
Gaussian  velocity  distributions  [47].  Numerical  solutions  of  the  10-  and  35-moment 
equations  have  been  obtained  for  one-dimensional  shocks  [14].  However,  the  10- 
moment  equations  do  not  include  heat-transfer  effects,  while  solutions  to  higher 
order  moment  equations,  including  the  13-,  14-,  20-,  21-  and  35-moment  equations, 
result  in  embedded  discontinuities  in  the  shock  structure  for  inflow  Mach  numbers 
higher  than  about  5  or  less  [14,  41,  74,  101].  Thus,  while  these  equations  are  valid  for 
higher  Knudsen  number  flows,  their  practical  utility  for  hypersonic  flows  is  limited. 

2.5  Simulation  Methods 

The  governing  equations  for  gas  flows  (the  Boltzmann  equation  and  the  Navier- 
Stokes  equations)  have  now  been  reviewed  in  general.  This  section  will  describe  two 
methods  used  to  simulate  gas  flows;  the  direct  simulation  Monte  Carlo  method  and 
Computational  Fluid  Dynamics. 

2.5.1  The  Direct  Simulation  Monte  Carlo  Method 

Although  the  Boltzmann  equation  is  valid  for  all  flow  regimes,  it  is  impossible  to 
solve  analytically  (except  for  extremely  simple  flows — although  analytical  solutions 
do  exist  for  collisionless  flows).  Numerically  solving  the  equation  quickly  becomes 
intractable  due  to  its  multi- dimensional  nature  (one  in  time,  three  in  physical  space 
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and  three  in  velocity  space)  and  the  complexity  of  the  collision  integral  term.  The 
direct  simulation  Monte  Carlo  (DSMC)  method  [4]  is  a  way  to  emulate  the  physical 
processes  modeled  by  the  Boltzmann  equation.  The  DSMC  method,  similar  to  other 
Monte  Carlo  schemes,  is  a  statistical  approach.  Instead  of  simulating  each  individ¬ 
ual  particle  in  a  gas  flow,  a  representative  sample  of  particles  is  followed  through  the 
flow.  Each  particle  has  a  specific  position,  velocity  and  internal  energy  (including  ro¬ 
tational  and  vibrational).  The  intermolecular  collisions  are  treated  on  a  probabilistic 
rather  than  a  deterministic  basis  and  assume  “molecular  chaos”  of  dilute  gas  flows. 
The  resulting  process  has  been  shown  to  converge  to  a  solution  of  the  Boltzmann 
equation  [95]. 

A  DSMC  implementation  can  be  briefly  described  as  follows.  A  physical  flow  do¬ 
main  with  appropriate  boundaries  is  described.  The  computational  domain  is  divided 
into  cells  used  for  selecting  collision  partners  and  over  which  the  particle  properties 
are  averaged  to  obtain  macroscopic  properties.  The  physical  domain  is  initialized 
with  a  number  of  representative  computational  particles  with  an  initial  position  and 
velocity  (according  to  an  equilibrium  VDF).  The  simulation  then  proceeds,  stepping 
through  time.  At  each  time  step 

•  The  particles  are  moved  according  to  the  velocity  and  time  step  size. 

•  Boundary  conditions,  such  as  collisions  with  walls,  inflow  and  outflow,  are 
applied. 

•  Particle  collisions  are  computed  based  on  collision  probabilities  and  molecular 
models. 

•  Macroscopic  properties  are  evaluated  by  taking  the  averages  of  the  properties 
of  the  individual  particles. 
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This  procedure  implies  certain  assumptions  and  limitations.  First,  the  time  step 
must  be  small  enough  relative  to  the  mean  collision  time  such  that  the  particle 
movements  and  the  collision  routines  can  be  separated.  Typical  limits  require  a 
time  step  to  be  approximately  1/3  of  the  mean  collision  time.  Second,  the  collision 
partners  are  chosen  based  on  the  particles  in  each  cell.  That  is,  each  cell  should 
be  less  than  one  mean  free  path  in  size — collision  partners  can  then  be  randomly 
chosen  from  the  particles  in  the  each  cell  while  maintaining  physical  accuracy.  Third, 
each  cell  should  contain  sufficient  particles  such  that  the  macroscopic  averages  are 
statistically  meaningful — 20  particles  per  cell  is  generally  required. 

DSMC  is  an  attractive  way  to  simulate  complex,  nonequilibrium  flows.  It  has 
been  shown  to  converge  to  solutions  of  the  Boltzmann  equation  in  the  limit  of  an  in¬ 
finite  number  of  particles  [95].  Both  the  Boltzmann  equation  and  DSMC  are  based  on 
the  same  physical  reasoning,  and  both  require  models  to  describe  surface  and  inter- 
molecular  interactions.  Nevertheless,  it  is  easier  to  implement  models  that  have  been 
phenomenologically  derived  to  agree  with  physical  reality  into  DSMC,  rather  than 
into  the  mathematically  rigorous  Boltzmann  equation  [84],  However,  the  practical 
utility  of  DSMC  is  limited  due  to  the  computational  cost.  As  the  Knudsen  number 
of  a  flow  decreases,  the  number  of  cells  (and,  hence,  particles)  required  increases. 
DSMC  simulations  of  higher  density  flows  are  limited  based  on  the  computational 
resources  available.  Thus,  DSMC  is  appropriate  for  the  simulation  of  flows  with  all 
types  of  nonequilibrium  in  the  transitional  and  rarefied  regimes. 

The  Lattice  Boltzmann  Equation 

The  Lattice  Boltzmann  Equation  is  a  “hyper-stylized  version  of  the  Boltzmann 
equation  explicitly  designed  to  solve  fluid-dynamics  problems”  [87].  LBE  methods 


are  suitable  for  flow  conditions  where  considerable  nonequilibrium  effects  are  present, 
such  as  rarefied  flows.  However,  these  methods  are  not  yet  suitable  for  flows  with 
strong  compressibility  and  substantial  heat  transfer  effects  [87];  their  applicability 
for  hypersonic  flows  is  thus  limited. 

2.5.2  Computational  Fluid  Dynamics 

The  Navier-Stokes  equations  can  be  solved  analytically  for  simple  flows,  or,  as 
is  the  case  for  Computational  Fluid  Dynamics  (CFD)  applications,  solved  numeri¬ 
cally.  (Strictly  speaking,  CFD  techniques  can  be  applied  to  any  of  the  conservation 
equations  previously  mentioned,  such  as  the  Burnett  or  13-moment  equations;  for 
this  study,  however,  it  should  be  noted  that  the  term  “CFD”  refers  to  numerically 
solving  the  Navier-Stokes  equations.)  The  finite- volume  method  is  commonly  used 
today  [37,  38,  89]. 

A  two-dimensional,  finite  volume  method  that  considers  a  single  species  would 
solve  the  Navier-Stokes  equations  in  conservative  form  as 

9Q  9(E,  -  E„)  9(F,  -  F„) 

at-  dx  +  dy  '  1  '  ’ 

Some  degree  of  vibrational  thermal  nonequilibrium  can  be  modeled  using  an  ad¬ 
ditional  energy  equation  for  the  vibrational  modes  [68].  The  rotational  modes  are 
assumed  to  be  in  equilibrium  with  the  translational  modes  and  are  modeled  with  one 
temperature,  T ;  the  vibrational  modes  are  modeled  with  a  vibrational  temperature, 
Tvib.  In  this  way,  the  vector  of  conserved  variables,  Q,  and  the  source  vector,  S,  are 
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given  as 
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where  p  is  the  mass  density,  u  and  v  are  the  bulk  velocity  in  the  x  and  y  directions, 
e  is  the  total  energy  and  ev  is  the  vibrational  energy  per  unit  volume  of  the  gas  and 
wv  is  the  vibrational  energy  source  term  (modeled  using  the  Landau- Teller  model  for 
vibrational  relaxation  [93]).  The  inviscid  and  viscous  flux  vectors  in  the  ^-direction 
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The  inviscid  and  viscous  flux  vectors  in  the  ^/-direction  are  similar.  The  shear  stress 
tensor  and  heat  flux  vector  are  evaluated  assuming  Newton’s  law  of  friction  and 
Fourier’s  law,  as  given  in  equations  2.12  and  2.13,  along  with  Stokes’  hypothesis 
leading  to  //#  =  0.  The  system  of  equations  is  closed  using  the  perfect  gas  law, 


p  =  pRT. 

CFD  has  been  used  for  many  years  to  simulate  gas  flows  in  the  continuum  regime. 
With  the  use  of  additional  equations,  it  can  also  successfully  simulate  vibrational 
(and  chemical)  nonequilibrium.  Typical  CFD  methods  do  not  include  an  additional 
rotational  energy  equation,  nor  can  they  accurately  model  the  translational  none¬ 
quilibrium  present  as  the  flow  becomes  more  rarefied.  The  computational  cost  of  a 
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CFD  simulation  is  a  function  of  the  number  of  cells  required,  which  is  not  explicitly 
dependent  on  the  flow  regime  in  question.  Thus,  CFD  is  appropriate  for  modeling 
high-temperature,  high-speed  gas  flows  in  the  continuum  regime,  with  limited  success 
in  the  slip  and  transitional  regimes  (with  the  slip  boundary  conditions). 

2.6  Computational  Codes 

The  remaining  chapters  describe  numerical  simulations  of  hypersonic  gas  flows 
obtained  using  the  DSMC  and  CFD  methods.  The  two  codes  used  to  obtain  these 
results  are  MONACO  (DSMC)  and  LeMANS  (CFD),  which  are  briefly  described 
here. 

2.6.1  MONACO 

In  this  research,  DSMC  results  are  provided  from  the  MONACO  code  [23],  which 
is  a  general  2D/3D,  object-oriented,  cell-based,  parallel  implementation  of  the  DSMC 
method.  It  has  been  applied  to  many  hypersonic,  rarefied  flows  [88].  Several  differ¬ 
ent  molecular  models  can  be  used,  including  the  Variable  Hard  Sphere  (VHS)  and 
Variable  Soft  Sphere  collision  models  [4,  46].  It  also  includes  variable  vibrational 
[92]  and  rotational  [9]  energy  exchange  probability  models  to  model  the  temperature 
dependence  of  the  rotational  and  vibrational  collision  probabilities.  In  cases  where 
the  cell  sizes  are  larger  than  the  local  mean  free  path,  the  subcell  method  can  be 
used  to  select  particles  for  collisions  [4]  to  ensure  physical  accuracy. 

The  current  work  uses  the  VHS  model,  the  details  of  which  are  discussed  in 


Chapter  III. 
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2.6.2  LeMANS 

The  CFD  results  are  obtained  using  the  (Le)  Michigan  Aerothermodynamic 
Navier-Stokes  (LeMANS)  code,  developed  at  the  University  of  Michigan  for  the  sim¬ 
ulation  of  hypersonic  reacting  flow-fields  [75,  76,  77].  LeMANS  is  a  general  2D/3D, 
parallel,  unstructured,  finite-volume  CFD  code  capable  of  simulating  gases  in  ther¬ 
mal  and  chemical  nonequilibrium.  A  modified  Steger- Warming  flux  vector  splitting 
scheme  is  implemented  [54,  85] .  The  code  is  second-order  in  space  and  time,  using  a 
point-  or  line-implicit  time  integration  scheme  [75].  A  two-temperature  model  [68]  is 
used  to  account  for  the  nonequilibrium  between  the  vibrational  and  the  translational- 
rotational  modes,  with  the  energy  exchange  rates  modeled  using  the  Landau- Teller 
model  [93].  Standard  finite-rate  chemistry  models  for  reacting  air  using  5  species, 
7  species  and  11  species  allow  atmospheric  flow  simulations  in  conditions  causing 
dissociation  and  weak  ionization  [68].  Transport  properties  can  be  computed  using 
Wilke’s  [103]  and  Blottner’s  [7]  models.  Different  boundary  conditions,  including 
no-slip  and  slip  velocity  and  temperature  jump  are  enforced. 

Details  concerning  the  transport  properties,  the  Landau- Teller  model  and  the  slip 
boundary  conditions  used  in  the  current  work  are  discussed  in  Chapter  III. 

A  Note  on  Hybrid  Methods 

DSMC  is  an  attractive  and  applicable  method  for  simulating  gas  flows  at  large 
Knudsen  numbers,  although  it  becomes  computationally  expensive  in  the  continuum 
regime.  CFD  works  well  for  the  continuum  regime,  but  is  inaccurate  for  higher 
Knudsen  number  flows.  Many  times,  the  gas  flow  around  a  hypersonic  vehicle  will 
fall  into  a  range  of  flow  regimes;  for  example,  there  might  be  a  region  of  high  density 
gas  upstream  of  the  vehicle  where  the  continuum  hypothesis  holds,  but  downstream 
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in  the  wake  the  gas  is  sufficiently  rarefied  that  the  continuum  hypothesis  breaks 
down.  Additionally,  nonequilibrium  effects  tend  to  be  concentrated  in  regions  of 
large  gradients,  such  as  shock  waves  and  boundary  layers.  Although  DSMC  could 
be  used  to  simulate  the  entire  domain,  it  is  computationally  expensive  to  do  so. 

A  class  of  hybrid  DSMC-CFD  solvers  has  been  proposed  to  solve  this  problem 
[36,  71,  80,  94],  In  a  hybrid  method,  the  computational  domain  is  divided  into 
continuum  and  particle  regions,  with  DSMC  being  used  in  the  particle  regions  and 
CFD  in  the  continuum  regions.  Two  main  problems  are  associated  with  this  type  of 
hybrid  method.  First  is  the  problem  of  determining  where  the  regions  of  equilibrium 
and  nonequilibrium  are  located.  This  is  usually  accomplished  through  some  type  of 
continuum  breakdown  parameter  [11],  which  is  discussed  in  some  detail  in  Chapter 
111.  Second  is  the  manner  in  which  information  is  passed  between  the  CFD  and 
DSMC  regions.  Research  continues  in  this  area  [79,  82,  80,  98]. 


CHAPTER  III 


Comparing  Simulation  Results  from  the  DSMC 

and  CFD  Methods 


3.1  Introduction 

The  direct  simulation  Monte  Carlo  (DSMC)  method  and  the  Computation  Fluid 
Dynamics  (CFD)  method  differ  in  their  fundamental  approach  to  simulating  gas 
flows.  Although  these  approaches  are  based  on  different  assumptions  regarding  the 
amount  of  equilibrium  present  in  the  flow,  each  should  give  equivalent  results  in 
regimes  where  they  are  both  valid.  The  goal  of  the  current  study,  to  determine  the 
effects  of  nonequilibrium  on  the  simulation  results,  requires  that  the  physical  models 
encapsulated  within  the  framework  of  each  method  be  treated  in  a  similar  manner. 
Thus  it  is  assured  that  differences  in  the  results  are  due  to  the  fundamental  differences 
of  the  methods,  and  their  ability  to  capture  nonequilibrium  phenomena,  and  not  due 
to  differences  in  their  treatment  of  the  underlying  physical  models.  Such  physical 
models  include  transport  properties  (such  as  viscosity),  wall  boundary  conditions 
and  vibrational  relaxation.  This  chapter  will  discuss  the  relevant  physical  models 
and  the  manner  in  which  they  are  treated  in  each  simulation  method  such  that  they 
are  equivalent,  as  much  as  is  possible. 

As  was  discussed  previously,  the  DSMC  method  is  valid  in  regions  of  nonequil- 
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ibrium,  while  the  CFD  method  is  valid  only  for  small  departures  from  equilibrium. 
It  is  instructive  to  determine  when  and  where  the  nonequilibrium  present  in  a  par¬ 
ticular  gas  flow  exceeds  the  capabilities  of  CFD.  Therefore,  methods  appropriate  for 
quantifying  the  amount  of  nonequilibrium  present  in  a  gas  flow,  in  the  form  of  a 
continuum  breakdown  parameter,  are  also  discussed. 

Much  of  the  material  of  this  chapter  is  also  relevant  when  designing  a  hybrid 
code  using  both  DSMC  and  CFD  methods.  In  a  hybrid  code,  it  is  critical  that  the 
physical  models  be  treated  in  a  consistent  manner  across  the  computational  domain. 
It  is  also  necessary  to  determine  where  in  the  computational  domain  DSMC  should 
be  used  rather  than  CFD — the  continuum  breakdown  parameter  is  used  in  that  role 
in  many  hybrid  methods. 

3.2  Transport  Properties 

Even  a  small  departure  from  LTE  in  a  gas  flow  due  to  variations  of  macroscopic 
phenomena,  such  as  species  concentration,  velocity  and  temperature,  will  give  rise  to 
the  transfer  of  mass,  momentum  and  energy,  respectively.  The  molecular  processes 
that  govern  the  transfer  of  mass,  momentum  and  energy  result  in  the  macroscopic 
properties  of  diffusion,  viscosity  and  thermal  conductivity  [93].  As  mentioned  above, 
the  goal  of  comparing  DSMC  and  CFD  simulations  requires  equivalent  treatment  of 
these  transport  properties.  As  all  simulations  treated  herein  involve  only  one  gas 
species,  the  transport  of  mass  via  diffusion  is  not  present  and  will  not  be  treated 
further  here. 

3.2.1  Viscosity 

All  DSMC  simulations  presented  here  are  generated  using  the  variable  hard  sphere 
(VHS)  model  [4],  The  viscosity  can  be  calculated  using  the  VHS  model  parameters. 
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Table  3.1:  Variable  hard  sphere  (VHS)  model  parameters  for  argon  and  nitrogen 
used  in  the  computational  simulations. 


Species 

u 

Tref  [K] 

dref 

H 

Ar 

0.734 

1000 

3.595  x  10~10 

n2 

0.7 

290 

4.110  x  10~10 

Therefore,  for  the  purposes  of  this  study,  LeMANS  was  modified  to  use  a  viscosity 
directly  corresponding  to  the  VHS  model  according  to  Eqs.  3.1  and  3.2  [81]: 


=  (3.1) 

l5\/irmkTref 

^ ref  =  27rd2ref  (5  -  2u)  (7  -  2 u) 

where  /j  is  the  viscosity,  T  is  the  temperature,  to  is  the  VHS  temperature  exponent, 
m  is  the  mass  of  one  molecule  of  the  gas,  k  is  the  Boltzmann  constant  and  d  is 
the  molecular  diameter.  This  model  requires  only  that  a  reference  temperature, 
Tref,  reference  diameter,  dref  and  the  temperature  exponent,  u,  be  specified  for 
both  DSMC  and  CFD  and  the  viscosity  treatment  is  equivalent,  for  a  gas  at  or  near 
equilibrium.  The  VHS  parameters  used  in  the  current  study  are  summarized  in  Table 
3.1. 

3.2.2  Thermal  Conductivity 


The  thermal  conductivity  coefficient,  k,  is  related  to  the  viscosity,  /i,  as 


Pr  = 


Cp/Jj 

K 


(3.3) 


where  Pr  is  the  Prandtl  number  and  cp  is  the  specific  heat  at  constant  pressure.  From 
kinetic  theory  an  equivalent  form  of  the  Prandtl  number  can  be  obtained,  known  as 


Eucken’s  Relation  [93], 
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Pr  = 


4y 


97-5 

from  which  the  thermal  conductivity  coefficient  is  obtained  as 


(3.4) 


97  -  5 

k  =  - R 11  .  (3. 

47-4  v 

The  ratio  of  specihc  heats,  7,  is  related  to  the  number  of  degrees  of  freedom,  £, 
present  in  a  particular  gas  species  as 


7  = 


£  +  2 
£ 


(3.6) 


The  use  of  Eqs.  3.5  and  3.6  in  LeMANS  ensures  equivalent  treatment  of  the  thermal 
conductivity  coefficient  in  both  methods,  for  a  gas  at  or  near  equilibrium. 


3.3  Wall  Boundary  Conditions 


Appropriate  treatment  of  the  boundary  conditions,  especially  the  wall  boundary 
conditions,  is  critical  when  comparing  solutions  from  DSMC  and  CFD.  This  sec¬ 
tion  discusses  the  different  wall  boundary  conditions  implemented  in  MONACO  and 
LeMANS. 

3.3.1  Gas-Surface  Interactions 

When  a  gas  molecule  collides  with  a  wall,  it  is  reflected  at  a  specihc  angle  and 
velocity.  In  DSMC,  this  effect  is  modeled  as  if  the  particle  were  absorbed  by  the  wall 
and  then  re-emitted  from  the  wall,  and  results  in  two  general  types  of  interactions. 

If  the  particle  is  reflected  elastically,  with  a  normal  velocity  component  equal  to 
the  negative  of  the  incoming  normal  velocity  component,  and  with  an  unchanged 
tangential  velocity  component,  the  collision  is  said  to  be  a  specular  interaction  [66]. 
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This  boundary  condition  corresponds  to  an  infinitely  smooth,  adiabatic  wall.  The 
wall  shear  stress  and  heat  transfer  rate  in  such  a  case  are  zero. 

The  other  effect  is  when  the  particle  collides  with  the  wall,  and  is  then  re-emitted 
at  a  random  angle  and  velocity  corresponding  to  a  Maxwellian  velocity  distribution 
in  equilibrium  with  the  wall.  The  colliding  particle  is  thus  thermally  accommodated 
to  the  wall.  This  collision  is  said  to  be  a  diffuse  interaction.  This  boundary  condition 
corresponds  to  a  rough  wall  held  at  a  constant  temperature.  The  wall  shear  stress 
and  heat  transfer  rate  are  correspondingly  non-zero. 

Both  of  these  conditions  are  theoretical  in  nature;  in  practice,  particles  collide 
with  walls  in  both  specular  and  diffuse  interactions.  The  accommodation  coefficient 
represents  the  fraction  of  incoming  particles  that  are  reflected  diffusely.  The  remain¬ 
der  are  reflected  specularly.  The  accommodation  coefficient  is  typically  a  function 
of  the  wall  roughness,  composition  and  the  types  of  surface  chemical  reactions  that 
might  occur. 

For  the  purposes  of  this  study,  only  the  case  of  a  fully  diffuse  interaction  is 
considered;  that  is,  the  accommodation  coefficient  is  unity,  and  the  walls  are  assumed 
to  be  isothermal,  with  all  temperatures  (translational,  rotational  and  vibrational) 
being  set  to  the  constant  wall  temperature. 

3.3.2  Velocity  Slip  and  Temperature  Jump 

At  sufficiently  low  Knudsen  numbers  to  warrant  the  assumption  of  continuum 
flow,  the  velocity  and  temperature  of  the  gas  near  the  wall  are  in  equilibrium  with 
the  wall.  That  is,  the  no-slip  boundary  conditions  hold.  As  the  flow  becomes  more 
rarefied,  there  are  insufficient  collisions  near  the  wall  to  equilibrate  the  gas  molecules 
to  the  wall,  thus  invalidating  the  no-slip  conditions.  This  results  in  conditions  known 
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as  slip  conditions  where  the  velocity  and  temperature  of  the  gas  near  the  wall  are  not 
equal  to  the  wall  velocity  and  temperature.  The  slip  conditions  are  inherent  to  the 
Boltzmann  equation,  and  are  thus  handled  transparently  by  the  DSMC  method.  It  is 
necessary,  however,  to  extend  the  no-slip  boundary  conditions  typically  used  in  CFD 
to  include  the  slip-conditions  to  improve  its  accuracy  when  used  for  more  rarefied 
flows.  This  sections  discusses  several  different  types  of  slip  boundary  conditions 
implemented  in  the  LeMANS  code,  as  well  as  the  method  to  extract  slip  quantities 
from  MONACO  simulations. 

Slip  Boundary  Conditions  in  CFD 


The  slip  boundary  condition  was  first  derived  by  Maxwell  for  a  flat  plate.  The 
expression  he  derived  for  the  velocity  slip  is  shown  in  Eq.  3.7,  as  discussed  in  Ref. 
[48],  where  wall  coordinates  are  used  (n  is  normal  to  the  wall  and  x  is  parallel  to  the 
wall): 


u.  =  a\2-^ + 


a  J  '  dn  '  4  pT  dx  ^  ^ 

where  Us  is  the  velocity  slip  (assuming  a  stationary  wall),  A  is  a  constant  of  propor¬ 
tionality,  a  is  the  tangential  momentum  accommodation  coefficient,  A  is  the  mean 
free-path,  ux  is  the  velocity  in  the  x-direction,  p  is  the  viscosity,  p  is  the  mass  density 
and  T  is  the  temperature.  For  an  isothermal  wall,  the  temperature  gradient  can  be 
neglected,  giving  the  boundary  condition  in  its  simplest  form  as: 


tt  a  i  2  -  (A  .  du 

U.  =  A  [  -  A 


a 


dn 


(3,8) 


n= 0 


The  mean  free  path,  A,  is  calculated  from  typical  gas  flow  properties  as  [29]: 


A  =  =  /W  tt 


pc  p  V  2  RT 


(3.9) 
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where  /i  is  the  viscosity,  p  is  the  mass  density  and  c  is  the  mean  molecular  speed. 
The  boundary  condition  for  the  temperature  jump  is  similar  [30]: 


_  2  -  a  27  dT 
a  (y  +  l)  Pr  dn  n=0 


(3.10) 


where  Tw  is  the  wall  temperature,  T0  is  the  temperature  of  the  gas  at  the  wall  (and 
where  Tw  —  T0  is  the  temperature  jump),  a  is  the  thermal  accommodation  coefficient, 
7  is  the  specific  heat  ratio,  Pr  is  the  Prandtl  number,  and  T  is  the  temperature.  An 
equivalent  mean  free  path  can  be  adopted  as 


At  — 


(7  +  1)  pccv  (7  +  1)  pcv 


(3.11) 


to  give  a  simpler  form, 

T0-Tw  =  —  At^  •  (3.12) 

«  9n  n= 0 

The  simple  Maxwell  slip  conditions  are  implemented  in  LeMANS  as  given  by  Eqs. 
3.8  and  3.12. 


The  general  form  of  the  velocity  slip  has  remained  the  same  throughout  many 
years  of  research;  the  difficulty  is  in  determining  the  correct  value  for  the  constant, 
A,  and  the  accommodation  coefficients,  a  and  a.  In  effect,  any  change  in  A  simply 
changes  the  value  of  the  accommodation  coefficient,  which  depends  on  the  physical 
characteristics  of  the  wall  itself  [86].  Many  times  these  coefficients  are  determined 
empirically,  while  other  times  they  are  calculated  from  relations  based  on  kinetic 
theory  [83].  The  most  detailed  boundary  conditions  require  knowledge  a  priori  of 
a  good  nonequilibrium  solution  for  the  flow  considered.  The  current  investigation 
seeks  to  find  a  good  boundary  condition  that  does  not  require  a  priori  knowledge 
of  any  nonequilibrium  solution,  which  would  eliminate  any  advantage  in  using  CFD 
over  DSMC.  For  the  purposes  of  this  study,  it  is  sufficient  that  the  values  for  the 
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Figure  3.1:  Velocity  profiles  in  the  Knudsen  layer.  The  true  velocity  profile  is  com¬ 
pared  to  Navier-Stokes  profiles  with  a  true  velocity  slip  and  a  fictitious 
velocity  slip  (after  Ref.  [49]). 

accommodation  coefficient  be  consistent  between  DSMC  and  CFD;  therefore,  an 
accommodation  coefficient  of  unity  is  used. 

The  previous,  simple  slip  boundary  conditions  were  derived  for  small  Knudsen 
numbers.  Gokgen  [29]  showed  that,  for  large  Knudsen  numbers,  these  simple  slip 
boundary  conditions  converge  to  different  values  than  those  predicted  by  free  molec¬ 
ular  flow,  ffe  then  proposed  the  general  slip  boundary  conditions, 

Q 

(®A  ®u>)  2A a  —  .  (3.13) 

Cln  n=  0 

Here,  a  is  either  velocity  or  temperature,  and  Aa  is  given  by  either  Eq.  3.9  or  3.11. 
For  small  Knudsen  numbers,  this  expression  reduces  to  Eqs.  3.8  and  3.12  [30]. 

The  Navier-Stokes  (N-S)  equations  with  no-slip  boundary  conditions  fail  in  at 
least  two  different  areas  as  the  flow  becomes  more  rarefied.  The  first  area,  the  as¬ 
sumption  of  no-slip  flow,  is  corrected  through  the  use  of  the  slip  boundary  conditions. 
However,  the  N-S  equations  also  assume  that  the  shear  stress  varies  linearly  with  the 
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velocity  gradient.  Cercignani  showed  that  this  is  not  the  case,  and  he  computed  the 
velocity  profile  near  the  wall,  as  discussed  in  [49].  Figure  3.1  illustrates  the  velocity 
profile  near  the  wall  in  the  Knudsen  layer  (which  is  defined  as  being  on  the  order  of 
one  mean  free-path  from  the  wall).  The  velocity  profiles  corresponding  to  the  N-S 
assumption  of  a  shear  stress  that  varies  linearly  with  the  velocity,  with  two  different 
values  of  velocity  slip,  are  shown  in  comparison  to  the  actual  velocity  profile. 

Here  there  are  two  choices  one  can  make  concerning  the  boundary  conditions;  use 
the  actual  velocity  slip  in  which  case  the  flow  field  away  from  the  wall  will  not  be 
accurate,  or  use  a  fictitious  velocity  slip  (as  shown  in  Figure  3.1).  The  fictitious  slip 
will  not  accurately  predict  the  flow  near  the  wall  (and  hence  may  adversely  affect 
the  accuracy  of  drag  and  heat  transfer  rates)  but  will  be  more  accurate  in  the  region 
away  from  the  wall. 

To  overcome  this  limitation,  a  correction  to  the  velocity  gradient  was  evaluated 
by  Lockerby,  et  al.  [49].  They  proposed  using  a  wall-function  type  of  boundary 
condition  for  the  velocity  where  the  viscosity  in  the  Knudsen  layer  is  modified  as 


where  the  wall-function, 

was  derived  from  a  curve-fit  to  Cercignani’s  Knudsen  layer  velocity  profile.  In  Eq. 
3.15,  n  is  the  distance  from  the  wall  and  A  is  the  mean  free-path. 

This  modification  of  the  viscosity  in  the  Knudsen  layer  is  also  used  in  connection 
with  a  simple  slip  boundary  condition  as  given  by  Eq.  3.8  with  A  =  \J 2/n.  This 
approach  is  expected  to  allow  the  CFD  method  to  accurately  model  the  velocity 
profile  in  the  Knudsen  layer. 
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This  new  wall  function  has  only  been  evaluated  for  isothermal  flow  conditions  that 
are  usually  encountered  in  micro- flows  [49].  Although  there  will  be  some  change  to 
the  heat  transfer  rate  coefficient  (based  on  a  constant  Prandtl  number),  there  might 
be  other  changes  necessary  to  give  the  correct  temperature  profile  in  the  Knudsen 
layer  for  the  non-isothermal,  hypersonic  flows  considered  here.  Nevertheless,  it  is 
instructive  to  investigate  the  possible  improvement  this  wall  function  might  afford. 
This  slip  condition  is  referenced  herein  as  the  Lockerby  slip  condition. 

For  cases  involving  a  diatomic  species  that  require  a  vibrational  energy  equation, 
a  corresponding  vibrational  temperature  jump  is  calculated  in  the  same  manner  as 
the  translational/rotational  temperature  jump,  with  the  appropriate  value  for  k  and 
cv  used  when  computing  the  mean  free-path  according  to  Eq.  3.11. 

The  CFD  solutions  with  velocity  slip  and  temperature  jump  boundary  conditions 
assume  a  fully  diffuse  wall  (a  =  1)  that  is  also  thermally  accommodating  (a  =  1). 
These  assumptions  correspond  to  the  fully  diffuse  wall  used  in  the  DSMC  simulations. 

Slip  Quantities  in  DSMC 

While  the  DSMC  method  inherently  includes  the  slip  boundary  conditions,  it  is 
necessary  to  extract  this  information  from  the  simulation  results.  In  MONACO,  the 
velocity  slip  and  temperature  jump  are  calculated  based  on  the  particles  that  strike 
the  surface,  according  to  Eqs.  3.16  and  3.17,  which  are  the  relations  used  in  Bird’s 
DSMC  implementation,  DS2V[5,  6]. 


EO/IV,!) 


(3.16) 


Ttra,j 


1  E  ((™/\Un\)  (||U||))  -  E  (m/\Un\)Us 
3  R  Z(V\Un\) 


Ttra,W 


(3.17) 
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where  Us  is  the  velocity  slip,  Un  is  the  velocity  normal  to  the  wall,  Up  is  the  velocity 
parallel  to  the  wall,  ||U||  is  the  velocity  magnitude  (includes  all  three  components), 
Ttraj  is  the  translational  temperature  jump,  and  Ttra,w  is  the  wall  translational 
temperature  and  the  summations  are  taken  over  all  particles  that  strike  the  surface. 

These  relations  include  the  somewhat  counterintuitive  use  of  the  velocity  com¬ 
ponent  normal  to  the  wall.  This  is  due  to  the  necessity  of  using  wall  fluxes,  an 
inherently  two-dimensional  quantity,  to  derive  a  velocity  slip  and  temperature  jump, 
which  are  volume  quantities.  That  is,  only  the  properties  of  particles  that  hit  the 
wall  are  considered,  rather  than  those  of  all  the  particles  in  the  cell  adjacent  to  the 
wall.  The  probability  that  a  particle  in  a  cell  will  collide  with  the  wall  during  a 
particular  time  step  is  proportional  to  its  velocity  component  normal  to  the  wall. 
Thus,  the  velocity  slip  and  temperature  jump  quantities  would  be  biased  toward 
higher  velocity  particles  if  the  summations  were  not  weighted  by  the  magnitude  of 
the  normal  velocity  component. 

Corresponding  quantities  can  be  derived  for  the  rotational  and  vibrational  tem¬ 
perature  jumps  as 


Trot,j 


i  E(<wm.p 

R  £(i/IK.|) 


Trot,W 


(3.18) 


qr  _ 

vib,j 


0, 


In  (  +  1 

f-'vib 


—  T, 


vib.W  j 


Tyib 


E  (<WI«»I) 


(3.19) 


EVV..I 

where  erot  and  evn  are  the  rotational  and  vibrational  energies  associated  with  each 
particle  and  ©„  is  the  characteristic  vibrational  temperature. 


A  Note  on  Temperatures 


As  shown  in  Chapter  II,  LeMANS  uses  a  two-temperature  model  (translation¬ 
al/rotational  and  vibrational)  while  MONACO  tracks  all  three  energy  modes:  trans- 
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lational,  rotational  and  vibrational.  Therefore,  the  proper  way  to  compare  simulation 
results  for  the  temperatures  is  to  average  the  DSMC  translational  and  rotational  tem¬ 
peratures  according  to  the  number  of  degrees  of  freedom  of  each  mode  (£tro  =  3  and 
=  2)  as  shown  in  Eq.  3.20  [4],  This  averaged  temperature  is  then  equivalent  to 
the  translational/rotational  temperature  obtained  from  LeMANS. 


3  Ttra  +  2  Trot 
5 


(3.20) 


In  the  same  way,  the  rotational  and  translational  temperature  jumps  obtained  from 
MONACO  are  averaged  and  then  compared  with  the  temperature  jump  predicted 
by  LeMANS. 

The  vibrational  temperature  results  obtained  from  each  simulation  method  are 
compared  directly. 


3.4  Vibrational  Relaxation 


The  rate  of  energy  exchange  between  the  translational  and  vibrational  modes 
is  inversely  proportional  to  the  vibrational  relaxation  time,  tv.  The  Landau- Teller 
model  is  generally  used  to  approximate  the  vibrational  relaxation  time  [93].  Millikan 
and  White  [56]  correlated  experimental  data  with  the  Landau-Teller  model  to  obtain 
an  expression  for  the  vibrational  relaxation  time  as 

TLT  =  lexip{^  +  B)  (3-21) 

where  p  is  the  pressure,  T  is  the  temperature  and  A  and  B  are  constants  that  vary 
by  species. 

A  correction  proposed  by  Park  [68]  that  is  necessary  for  temperatures  typically 
encountered  at  hypersonic  speeds  is  given  as 

1 

rP=  — 


nca 


(3.22) 
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where  n  is  the  particle  density,  c  is  the  mean  molecular  speed  and  a  is  the  vibrational 
collision  cross  section. 

The  vibrational  relaxation  time  is  then  given  in  Eq.  3.23. 

tv  =  tit  +  tp  (3.23) 

The  Landau- Teller  model  is  a  continuum  model  and  can  be  directly  implemented 
in  a  CFD  code.  The  vibrational  relaxation  process  in  DSMC  is  governed  by  collision 
probabilities.  The  probability  that  a  collision  will  result  in  an  exchange  of  vibrational 
energy  is  given  as 

P  =  —  (3.24) 

tvv 

where  v  is  the  bi-molecular  collision  rate.  This  probability  is  then  implemented 
in  MONACO  using  discrete  collision  probabilities  [10,  92]  (with  various  corrections 
[53,  1]),  which,  when  integrated  over  all  collisions,  is  expected  to  correspond  to  the 
total  probability  as  given  in  Eq.  3.24. 

Proper  comparison  of  the  vibrational  relaxation  in  both  CFD  and  DSMC  then 
requires  the  use  of  equivalent  constants,  A  and  B ,  in  the  Millikan  and  White  correla¬ 
tion  for  vibrational  relaxation  time  as  well  as  equivalent  Park  correction  parameters. 
The  values  for  A  and  B  for  nitrogen  are  220  and  -24.80,  respectively,  for  a  pressure 
given  in  atm.  For  a  pressure  given  in  Pa,  the  value  for  B  is  -13.27.  LeMANS  uses 
the  Millikan  and  White  values  in  the  Landau- Teller  model,  although  with  a  different 
definition  of  B  that  corresponds  directly  to  B  —  -24.8.  The  current  work  uses  a 
constant  vibrational  collision  cross  section,  a  =  5.81  x  ICC21  m2,  for  both  the  CFD 
and  DSMC  simulations. 

An  additional  correction  to  the  vibrational  collision  probability  in  MONACO  is 
necessary  to  achieve  accurate  results.  A  recent  study  compared  the  theoretical  vi- 
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Table  3.2:  Vibrational  collision  probability,  P ,  in  MONACO  compared  to  theory  for 


N2-N2  collisions,  with  the  correction  factor  [24], 


Temperature  (K) 

MONACO 

Theory 

Correction  Factor 

5000 

1.57  x  10“4 

1.24  x  10“4 

0.79 

10000 

1.71  x  10“3 

2.44  x  10"3 

1.43 

15000 

4.20  x  10“3 

7.68  x  10"3 

1.83 

20000 

6.76  x  10“3 

1.23  x  10“2 

1.82 

25000 

9.00  x  10“3 

1.53  x  10“2 

1.70 

30000 

1.08  x  10~2 

1.72  x  10“2 

1.59 

35000 

1.24  x  10~2 

1.85  x  10“2 

1.49 

40000 

1.38  x  10“2 

1.95  x  10“2 

1.41 

45000 

1.51  x  10“2 

2.03  x  10“2 

1.34 

50000 

1.61  x  10"2 

2.10  x  10“2 

1.30 

■e -  MONACO 


Temperature  [K] 
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Figure  3.2:  Vibrational  collision  probability  for  N2-N2  collisions  in  MONACO  com¬ 
pared  to  theory,  with  the  correction  factor. 
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brational  collision  probability  (Eq.  3.24)  with  that  given  by  MONACO  [24],  This 
study  showed  that,  while  remaining  within  a  factor  of  2,  the  overall  probability  seen 
in  MONACO  differed  from  the  theoretical  probability  as  shown  in  Table  3.2  and 
Figure  3.2  for  N2-N2  collisions.  The  correction  factor  is  the  ratio  of  the  theoreti¬ 
cal  vibrational  collision  probability  to  that  of  MONACO.  The  discrepancies  between 
MONACO  and  the  theoretical  vibrational  probability  are  most  likely  due  to  the 
method  of  analytical  integration  used  in  deriving  the  discrete  collision  probabili¬ 
ties  used  in  the  DSMC  model  [10].  It  should  be  emphasized  that  this  correction 
is  obtained  by  comparing  the  vibrational  collision  probability  in  MONACO  to  the 
theoretical  vibrational  collision  probability,  and  not  to  any  CFD  results. 

In  the  present  work,  the  vibrational  collision  probability  is  multiplied  by  a  cor¬ 
rection  factor  (as  shown  in  the  fourth  column  of  Table  3.2)  to  better  correspond  to 
theory.  The  correction  factor  is  chosen  based  on  the  maximum  translational  tem¬ 
perature  in  each  case;  for  example,  in  a  Mach  10  flow  the  maximum  temperature  is 
approximately  5,000  K,  so  the  vibrational  collision  probability  is  multiplied  by  0.79, 
while  for  a  Mach  25  case  the  maximum  temperature  is  approximately  25,000  and  the 
vibrational  collision  probability  is  multiplied  by  1.70. 

3.5  Continuum  Breakdown/Nonequilibrium  Onset 

The  areas  of  the  flow  where  the  continuum  hypothesis  breaks  down  (or  equiv¬ 
alently,  where  the  flow  is  no  longer  in  local  thermodynamic  equilibrium),  can  be 
quantified  by  the  use  of  a  continuum  breakdown  parameter.  Several  breakdown  pa¬ 
rameters  have  been  proposed  in  the  literature.  Bird  [3]  proposed  the  parameter,  P, 
for  jet  expansions: 

P=LL  dJL  =mM  1  ^ 

pv  ds  V  8  P  ds 


(3.25) 
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where  U  is  the  local  velocity,  v  is  the  collision  frequency,  M  is  the  local  Mach  number, 
7  is  the  ratio  of  specific  heats,  and  s  is  the  distance  along  a  streamline.  Continuum 
breakdown  is  assumed  when  P  is  greater  than  0.05. 

Boyd,  et  al.  [11]  carried  out  extensive  numerical  investigation  of  one-dimensional 
normal  shock  waves  and  two-dimensional  bow  shocks  comparing  DSMC  and  CFD 
results  to  determine  an  appropriate  breakdown  parameter.  They  proposed  the  use 
of  the  gradient-length  local  Knudsen  number, 


Kugll  =  g 


dQ 

~di 


(3.26) 


where  A  is  the  mean  free-path,  Q  is  some  quantity  of  interest  such  as  density,  pressure, 
temperature  or  velocity  magnitude,  and  the  derivative  is  taken  in  the  direction  of 
the  maximum  gradient. 

Other  breakdown  parameters  that  have  been  proposed  include  Tiwari’s  crite¬ 
rion  [90], 


Id’ll  = 


2^+1-\\r\ 
pRT  \  5  RT  2  11  1 


(3.27) 


where  p  is  the  mass  density,  R  is  the  gas  constant,  T  is  the  temperature,  ||r||£  is 
the  Euclidean  norm  of  the  stress  tensor,  r,  and  q  is  the  heat  flux  vector;  and  the 
parameter  B  [26], 


B  =  max  { | T*j |  ,  \q*\} 


(3.28) 


where  r,*-  and  q*  are  the  normalized  stress  tensor  and  heat  flux  vector,  respectively. 
More  recent  work  [15]  has  proposed  the  use  of  entropy  generation  rates  to  determine 
areas  of  nonequilibrium. 

Further  studies  [98]  have  shown  that  Bird’s  parameter  (based  on  density  and 
other  properties  such  as  temperature  and  velocity  magnitude)  is  inadequate  in  the 
stagnation  region  of  hypersonic  compressible  flows.  Of  the  alternatives,  only  the 
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gradient-length  local  Knudsen  number,  Kiigll,  has  been  extensively  tested  for  hy¬ 
personic  compressible  flows.  Therefore  this  study  will  use  Kiigll,  computed  using 
Eq.  3.26,  as  the  breakdown  parameter. 

When  calculating  KiiGll  based  on  velocity  magnitude,  the  gradient  is  normalized 
by  the  maximum  of  the  local  velocity  magnitude  and  the  local  speed  of  sound,  ft 
is  generally  assumed  that  continuum  breakdown  occurs  whenever  Kiiqll  is  greater 
than  0.05. 


CHAPTER  IV 


Hypersonic  Flow  about  a  Cylinder 

4.1  Introduction 

Hypersonic  vehicles  can  be  broadly  classified  as  either  blnnt-body  or  sharp- 
leading-edge  vehicles.  The  flows  around  each  of  these  two  types  of  vehicles  are 
significantly  different  and  highlight  different  physical  phenomena.  This  chapter  con¬ 
siders  a  hypersonic  flow  about  a  blunt-body;  in  this  case  a  two-dimensional,  12-inch 
diameter  cylinder  as  shown  in  Figure  4.1.  Results  from  computational  simulations 
are  discussed. 

First,  the  case  of  a  hypersonic  flow  of  argon,  a  monatomic  gas,  is  considered.  The 
monatomic  nature  of  argon  eliminates  the  possibility  of  thermal  nonequilibrium  clue 
to  the  nonexistence  of  internal  energy  modes.  The  results  from  this  part  of  the  study 
will  help  to  understand  the  effects  of  translational  nonequilibrium  on  the  flow  field 
and  surface  property  predictions,  and  will  give  a  baseline  to  which  further  results 
can  be  compared.  Then  the  case  of  a  hypersonic  flow  of  nitrogen,  using  the  same 
physical  geometry,  is  considered.  The  use  of  nitrogen,  a  diatomic  gas,  allows  the 
investigation  of  the  effects  of  thermal  nonequilibrium  in  addition  to  the  translational 
nonequilibrium  present  in  the  argon  flow.  It  should  be  noted  that  only  thermal 
nonequilibrium  effects  are  considered;  although  the  temperatures  present,  especially 
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Figure  4.1:  2D  cylinder  geometry  definition.  <f>  is  the  angle  in  degrees  from  the 
stagnation  point. 

at  the  higher  velocities,  are  high  enough  that  significant  amounts  of  dissociation  and 
ionization  would  be  expected,  no  chemistry  effects  are  considered. 

The  inflow  boundary  conditions  are  Mach  10  and  Mach  25,  and  the  free  stream 
density  of  the  flow  is  varied  such  that  several  different  regimes  are  considered,  from 
the  continuum  through  the  transitional  to  the  rarefied  regime,  as  shown  in  Table 
4.1.  Here,  the  global  Knudsen  number  is  calculated  based  on  the  cylinder  diameter, 
using  the  hard-sphere  model  for  the  mean  free-path,  A, 

A  =  (4'1} 

where  d  is  the  molecular  diameter  and  n  is  the  particle  density.  Note  that  in  all  cases 
the  flow  remains  in  the  laminar  regime,  as  shown  by  the  small  Reynolds  number 
values. 

Other  relevant  in-flow  boundary  conditions  and  the  constant  wall  temperature 
are  shown  in  Table  4.2. 

Surface  and  flow  held  properties  for  this  how  are  presented  from  two  different 
computational  approaches.  First,  CFD  results  are  obtained  through  solution  of  the 
Navier-Stokes  equations,  using  the  LeMANS  code.  Different  boundary  conditions, 
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Table  4.1:  Flow  regimes  considered. 


w 

n  [particles/m3] 

Poo  (Ar)  [kg/m3] 

Poo  (N2)  [kg/m3] 

Re  ** 

0.002 

0.01 

0.05 

0.25 

2.124  x  1021 
4.247  x  1020 
8.494  x  1019 
1.699  x  1019 

1.408  x  10"4 
2.818  x  10“5 
5.636  x  10“6 
1.127  x  10“6 

9.872  x  10“5 
1.974  x  10“5 
3.949  x  10“6 
7.897  x  10“7 

20,000 

4,000 

800 

160 

Based  on  hard-sphere  mean  free  path. 
Based  on  cylinder  diameter  and  Mach  25. 


Table  4.2:  Boundary  conditions. 


Mach 

Uoo  (Ar)  [m/s] 

I/oo  (N2)  [m/s] 

Too  [K] 

Twall  [K] 

10 

2624 

2883 

200 

500 

25 

6585 

7208 

200 

1500 

including  no-slip  and  slip  velocity  and  temperature  jump  are  enforced,  as  discussed 
in  Chapter  III.  Recall  that  the  boundary  conditions  implemented  are  the  no-slip 
conditions,  the  simple  Maxwell  slip  conditions,  the  Gokgen  slip  conditions  [29]  and 
the  Knudsen-layer  correction  to  the  standard  slip  conditions  proposed  by  Lockerby, 
et  al.  [49]. 

In  each  case,  a  grid  independence  study  is  conducted  to  determine  the  final  mesh 
resolution  used.  As  the  current  study  is  most  interested  in  the  surface  properties, 
the  surface  property  profiles  and  integrated  quantities  are  used  to  define  a  mesh- 
indendent  solution.  Successive  solutions  are  generated  as  the  number  of  nodes  in 
the  wall-normal  direction  in  an  area  near  the  wall  containing  the  boundary  layer 
is  doubled,  and  the  number  of  nodes  in  the  wall-parallel  direction  is  doubled.  A 
mesh-independent  solution  is  defined  as  one  for  which  the  total  drag  and  peak  heat 
transfer  rate  change  by  1%  or  less,  and  there  is  no  discernable  change  in  the  overall 
surface  property  profiles  of  pressure,  shear  stress  and  heat  flux. 
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Second,  DSMC  results  are  provided  from  the  MONACO  code  for  the  same  flow 
conditions.  In  general,  the  mesh  used  for  the  final  solution  for  each  case  is  adapted 
from  previous  solutions  such  that  each  cell  size  is  on  the  order  of  a  mean  free  path  or 
smaller.  The  exceptions  are  the  Kn  =  0.002  cases,  where  the  cell  size  is  approximately 
four  times  the  mean  free  path,  and  the  Kn  =  0.01  cases,  where  the  cell  sizes  near  the 
surface  are  on  the  order  of  two  mean  free  paths  in  size.  In  these  cases  the  subcell 
method  is  used  to  select  particles  for  collisions  to  ensure  physical  accuracy  [4], 

Example  meshes  for  both  DSMC  and  CFD  for  the  different  flow  regimes  are 
shown  in  Figure  4.2. 

It  should  be  noted  here  that  the  DSMC  solutions  are  assumed  to  be  the  more 
correct  solutions,  as  discussed  in  Chapter  II. 

In  the  results  that  follow,  the  surface  properties  are  presented  in  terms  of  non- 
dimensional  coefficients, 


n  _  P-Poo 
P  poo  Ul 

CF=T  T 


Ch  —  i 


2  r°0 OO 

q 


(4.2) 

(4.3) 

(4.4) 


\PooUl 

where  p  is  the  pressure,  r  is  the  shear  stress,  q  is  the  heat  transfer  rate,  p^  is  the  free 
stream  pressure,  p ^  is  the  free  stream  density  and  is  the  free  stream  velocity. 
The  surface  properties  in  each  case  are  plotted  as  a  function  of  the  angle  around  the 
cylinder,  with  the  stagnation  point  being  located  at  an  angle  of  zero  (see  Figure  4.1). 

Along  with  the  surface  properties,  the  maximum  value  of  Kiigll  at  the  surface 
(based  on  the  CFD  solution)  is  plotted  in  each  case.  Note  that  the  re-axis  crosses  the 
right  y-axis  at  KncLL.crit  =  0.05,  thus  only  continuum  breakdown  parameter  values 


above  the  critical  value  are  shown  in  the  figures. 


(a)  Kn  =  0.002 


Figure  4.2:  Example  meshes  for  both 
Mach  10,  nitrogen  meshes 


(d)  Kn  =  0.25 
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(a)  Mach  10 


(b)  Mach  25 


Figure  4.3:  Percentage  of  total  drag  due  to  pressure  and  friction  for  a  flow  of  argon 
about  a  cylinder. 


4.2  Argon 


The  total  drag  and  peak  heat  transfer  rates  predicted  by  both  computational 
methods  for  a  hypersonic  flow  of  argon  around  the  cylinder  are  summarized  in  Tables 
4.3  and  4.4.  The  differences  between  CFD  and  DSMC  are  also  shown  graphically  in 
Figures  4.4  and  4.5.  Here,  the  peak  heat  transfer  rate  is  obtained  by  averaging  over 
the  surface  area  within  one  degree  of  the  stagnation  point.  For  CFD,  these  quantities 
are  calculated  for  each  of  the  different  boundary  conditions  implemented. 

Figure  4.3  illustrates  the  percentage  of  total  drag  due  to  pressure  and  friction 
forces,  for  both  DSMC  and  CFD.  It  is  significant  to  note  that  as  the  Knudsen  number 
increases,  the  percentage  of  total  drag  due  to  friction  increases  from  about  5%  at 
Kn  =  0.002  to  about  20%  at  Kn  =  0.25.  Also  note  that  the  vast  majority  of  the 
drag  is  due  to  pressure  forces.  It  is  also  significant  to  note,  as  shown  in  Figure  4.4, 
that  the  difference  in  predicted  total  drag  between  CFD  and  DSMC  is  due  mostly 
to  the  differences  predicted  in  the  friction  forces.  Also,  the  large  disagreement  with 
the  no-slip  condition  results  is  due  to  the  increase  in  predicted  friction  drag. 


Difference  from  DSMC 
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Table  4.3:  Total  drag  a  for  flow  of  argon  about  a  cylinder. 

Mach  10 


Drag/Length  [N/m]  (%  Difference) 

Kn^ 

DSMC 

CFD 

no-slip 

Maxwell 

Gokgen 

Locker  by 

0.002 

187.6 

187.5  (-0.1%) 

187.4  (-0.1%) 

187.4  (-0.1%) 

187.6  (0.0%) 

0.01 

40.02 

40.30  (0.7%) 

40.17  (0.4%) 

40.11  (0.2%) 

40.20  (0.5%) 

0.05 

8.866 

9.358  (5.6%) 

9.082  (2.4%) 

8.921  (0.6%) 

9.148  (3.2%) 

0.25 

2.092 

2.579  (23.3%) 

2.296  (9.7%) 

2.067  (-1.2%) 

2.363  (12.9%) 

Mach  25 


Krioo 

Drag/Length  [N/m]  (%  Difference) 

DSMC 

CFD 

no-slip 

Maxwell 

Gokgen 

Lockerby 

0.002 

1171 

1176  (0.4%) 

1176  (0.4%) 

1176  (0.4%) 

1177  (0.5%) 

0.01 

250.8 

255.3  (1.8%) 

254.5  (1.5%) 

254.3  (1.4%) 

254.8  (1.6%) 

0.05 

56.92 

60.82  (6.9%) 

58.89  (3.5%) 

57.63  (1.3%) 

59.38  (4.3%) 

0.25 

13.19 

17.66  (34.0%) 

15.64  (18.6%) 

13.94  (5.7%) 

16.18  (22.7%) 

40% 
35%  h 


<j  30%  h  ■  ■ 


V) 

Q 


25%  - 
20%  - 
15%  - 
10% 
5%  h 
0% 
-5% 


i  i  i - 1  i  i 


HI  Friction 
Pressure 


J 


i  I  I _ I  I  I 
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0.002  0.01  0.05  0.25 

Knudsen  Number 


0.002  0.01  0.05  0.25 

Knudsen  Number 


(a)  Mach  10 


(b)  Mach  25 


Figure  4.4:  Total  drag  difference  from  DSMC  for  a  flow  of  argon  about  a  cylinder. 

Note  that  there  is  a  negative  contribution  to  total  drag  due  to  friction 
for  the  Mach  25,  Kn  =  0.25  case,  bringing  the  total  drag  difference  down 
to  5.7%. 


Difference  from  DSMC 
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Table  4.4:  Peak  heat  transfer  rate  for  a  flow  of  argon  about  a  cylinder. 

Mach  10 


Peak  Heating  [kW / m2] 

(%  Difference) 

Klloo 

DSMC 

CFD 

no-slip 

Maxwell 

Gokgen 

Locker  by 

0.002 

89.80 

89.84  (0.1%) 

89.14  (-0.7%) 

89.00  (-0.9%) 

88.54  (-1.4%) 

0.01 

39.13 

39.76  (1.6%) 

39.06  (-0.2%) 

38.77  (-0.9%) 

39.23  (0.3%) 

0.05 

15.85 

18.02  (13.7%) 

17.20  (8.5%) 

16.63  (5.0%) 

17.43  (10.0%) 

0.25 

5.926 

7.828  (32.1%) 

7.040  (18.8%) 

6.011  (1.4%) 

7.296  (23.1%) 

Mach  25 


Peak  Heating  [kW / m2] 

(%  Difference) 

Klloo 

DSMC 

CFD 

no-slip 

Maxwell 

Gokgen 

Locker  by 

0.002 

1746 

1763  (0.9%) 

1750  (0.2%) 

1746  (-0.1%) 

1739  (-0.5%) 

0.01 

749.6 

791.0  (5.5%) 

771.1  (2.9%) 

762.1  (1.7%) 

774.8  (3.4%) 

0.05 

304.5 

357.9  (17.1%) 

339.5  (11.5%) 

321.2  (5.5%) 

345.0  (13.3%) 

0.25 

108.2 

148.2  (37.0%) 

133.4  (23.2%) 

106.5  (-1.6%) 

138.6  (28.0%) 

Knudsen  Number  Knudsen  Number 

(a)  Mach  10  (b)  Mach  25 

Figure  4.5:  Peak  heat  transfer  rate  difference  from  DSMC  for  a  flow  of  argon  about 
a  cylinder. 
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It  is  clear  that  the  slip  boundary  conditions  improve  the  agreement  between  the 
two  simulation  methods.  While  the  difference  in  total  drag  is  as  high  as  23.3%  and 
34.0%  for  the  no-slip  condition,  the  Gokgen  slip  boundary  condition  shows  the  best 
agreement,  with  the  difference  exceeding  5%  only  for  the  Mach  25,  Kn  =  0.25  case. 
Similar  trends  are  noted  for  the  peak  heat  transfer  rate,  with  a  37%  difference  for 
the  no-slip  CFD  solution,  while  the  Gokgen  slip  solution  remains  near  2%  or  less  for 
most  of  the  cases  and  is  5%  or  greater  only  for  the  Kn  =  0.05  cases. 

There  is  also  a  slight  improvement  in  drag  and  heat  flux  agreement  between  the 
Maxwell  and  Lockerby  boundary  conditions,  where  the  only  difference  in  the  two 
boundary  conditions  is  the  inclusion  in  the  Lockerby  slip  condition  of  the  viscosity 
correction  function  within  the  Knndsen  layer,  and  the  nse  of  a  more  correct  value 
for  the  constant  of  proportionality,  A,  in  Eq.  3.8.  Surprisingly,  it  is  the  Maxwell 
boundary  condition  that  shows  better  agreement  with  DSMC  than  the  Lockerby 
boundary  condition.  This  improvement  is  not  as  great  as  that  achieved  with  the 
Gokgen  boundary  condition,  but  there  is  value  in  using  this  simple  slip  boundary 
condition. 

Although  there  is  very  good  agreement  between  the  DSMC  and  Gokgen  CFD  re¬ 
sults  for  the  surface  properties,  there  are  still  significant  differences  in  flow  properties 
(such  as  shock  structure),  as  shown  below. 

Due  to  its  better  agreement  with  DSMC,  all  CFD  results  shown  below  are  taken 
from  the  Gokgen  CFD  cases  unless  otherwise  specified. 

4.2.1  Continuum  Breakdown 

The  breakdown  parameter  is  calculated  using  both  the  CFD  and  the  DSMC 
solutions  according  to  Equation  3.26,  with  the  derivative  being  taken  in  the  direction 
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of  the  steepest  gradient.  For  the  case  of  a  cylinder  in  a  hypersonic  flow  of  a  simple  gas, 
the  only  causes  of  breakdown  to  the  continuum  hypothesis  expected  are  in  regions 
of  high  gradients  (such  as  the  shock  and  boundary  layer)  and  regions  of  rarefaction 
(such  as  in  the  wake).  The  amount  of  continuum  breakdown  is  also  expected  to 
increase  as  the  gas  flow  becomes  more  rarefied. 

These  trends  are  confirmed  by  the  results  that  are  shown  in  Figures  4.6  and 
4.7,  which  illustrate  the  amount  of  continuum  breakdown  in  the  flow  as  it  becomes 
more  rarefied,  for  the  Mach  10  and  Mach  25  cases,  respectively.  The  maximum 
gradient  length  local  Knudsen  number  is  computed  from  the  DSMC  (top)  and  CFD 
(bottom)  solutions.  Light  gray  regions  correspond  to  KnGLL  <  0.05,  dark  gray 
regions  correspond  to  Kiiqll  <  0.10  and  black  regions  correspond  to  Kiiqll  >  0.10. 

In  general,  the  flow  experiences  continuum  breakdown  in  three  areas;  across  the 
bow  shock,  in  the  boundary  layer  and  in  the  wake  region.  The  flow  in  the  shock  and 
boundary  layer  experiences  very  steep  gradients  in  flow  properties,  while  the  wake 
region  is  more  rarefied,  thus  leading  to  the  breakdown  of  the  continuum  hypothesis. 

At  a  Knudsen  number  of  0.002,  the  flow  is  within  the  continuum,  no-slip  regime. 
Nevertheless,  there  is  still  evidence  of  continuum  breakdown  in  the  shock,  along  the 
cylinder  surface  in  the  boundary  layer  and  in  the  wake.  Interestingly,  DSMC  predicts 
a  larger  degree  of  breakdown  than  does  CFD.  A  Knudsen  number  of  0.01  is  considered 
to  be  near  the  limit  of  the  continuum,  no-slip  regime,  with  increased  evidence  of 
continuum  breakdown.  At  a  Knudsen  number  of  0.05,  the  flow  is  well  within  the  slip 
regime  with  an  associated  increase  in  the  regions  of  continuum  breakdown.  For  a 
Knudsen  number  of  0.25,  the  flow  would  be  considered  outside  of  the  slip  regime  and 
into  the  transition  regime.  Here,  even  the  addition  of  slip  boundary  conditions  is  not 
expected  to  help  the  continuum  CFD  method’s  predictive  capabilities  very  much. 
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Figure  4.6:  Kiigll  field  for  a  Mach  10  flow  of  argon  about  a  cylinder.  Light  gray 
regions  correspond  to  Kiigll  <  0.05,  dark  gray  regions  correspond  to 
Kiigll  <  0.10  and  black  regions  correspond  to  Kiiqll  >  0.10 


Figure  4.7:  Kiigll  field  for  a  Mach  25  flow  of  argon  about  a  cylinder.  Light  gray 
regions  correspond  to  Kiigll  <  0.05,  dark  gray  regions  correspond  to 
Kiigll  <  0.10  and  black  regions  correspond  to  Kiiqll  >  0.10 
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Indeed,  the  plots  of  the  breakdown  parameter  indicate  that  there  are  significant 
nonequilibrium  effects  across  almost  all  of  the  flow  domain. 

Similar  trends  to  those  seen  at  Mach  10  are  also  present  at  Mach  25;  the  main 
difference  being  that  continuum  breakdown  increases  more  rapidly  as  the  flow  be¬ 
comes  more  rarefied.  This  is  expected  as  the  gradients  in  flow  properties  are  expected 
to  be  more  severe;  as  the  flow  velocity  increases,  the  residence  time  of  a  fluid  par¬ 
ticle  decreases,  leaving  less  time  for  the  properties  to  equilibrate  to  the  local  flow 
conditions. 

Figures  4.6  and  4.7  can  also  give  an  idea  of  the  relative  efficiency  of  a  hybrid 
CFD/DSMC  method  using  Kiiqll  >  0.05  as  an  indicator  of  regions  where  the  flow 
can  no  longer  be  accurately  predicted  by  CFD.  With  such  a  design,  the  portions  of 
the  computational  domain  represented  by  the  dark  gray  and  black  regions  should 
be  computed  using  DSMC.  Initially,  it  would  be  necessary  to  base  the  interfaces 
according  to  the  value  of  the  continuum  breakdown  parameter  predicted  by  CFD 
(top  of  sub-figures).  As  the  solution  progresses,  the  final  interface  locations  should 
be  determined  by  the  more  accurate,  hybrid  solution,  which  would  correspond  to  the 
DSMC  solution  (shown  at  the  bottom  of  each  sub- figure).  Hence,  the  final  interfaces 
would  be  determined  by  the  amount  of  continuum  flow  predicted  by  the  DSMC 
solutions. 

With  this  in  mind,  it  can  be  seen  that  there  would  be  very  small  portions  of  the 
computational  domain  requiring  DSMC  for  the  Kn  =  0.002  and  Kn  =  0.01  cases, 
while  for  the  Kn  =  0.05  case  only  small  portions  remain  in  the  CFD  domain.  The 
amount  of  continuum  breakdown  present  in  the  Kn  =  0.25  case  implies  that  this 
would  not  be  a  likely  candidate  for  a  hybrid  method;  nearly  the  entire  domain  would 
require  DSMC. 
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4.2.2  Flow  Field  Properties 

The  density  ratio  fields,  where  the  density  is  normalized  by  the  freestream  density, 
can  be  seen  in  Figures  4.8  and  4.9.  Note  that  the  maximum  density  ratio  behind 
the  shock  is  about  4,  agreeing  with  the  theoretical  limit  for  a  monatomic  gas  in 
equilibrium  and  the  results  from  the  normal  shock  relations,  p/po  =  3.87  (Mach  10) 
and  p/p0  =  3.97  (Mach  25).  Overall  agreement  between  CFD  and  DSMC  is  good, 
with  some  small  differences  in  the  shock  structure  at  the  higher  Knudsen  number 
conditions.  Also  note  that  as  the  freestream  density  decreases,  the  shock  becomes 
much  more  diffuse. 

The  temperature  fields  predicted  by  both  CFD  and  DSMC  can  be  seen  in  Figures 
4.10  and  4.11  for  the  Mach  10  and  Mach  25  cases,  respectively.  The  Kn  =  0.002  tem¬ 
perature  field  exhibits  the  typical  flow  features  found  in  a  blunt  body  flow;  namely, 
a  fairly  thin  shockwave  standing  off  from  the  body,  a  high  temperature  region  fol¬ 
lowing  the  shock  (as  the  translational  energy  is  converted  into  thermal  energy  with 
the  decrease  in  velocity  due  to  shock  compression),  a  thermal  boundary  layer  that 
gradually  thickens  around  the  cylinder,  and  a  wake  of  lower  temperature  as  the  flow 
expands  around  the  top  of  the  cylinder.  The  maximum  temperature  behind  the 
shock  agrees  with  that  from  the  normal  shock  relations,  Tmax  =  6, 460  K  (Mach  10) 
and  Tmax  =  39,400  K  (Mach  25).  The  two  simulation  methods  agree  very  well  for 
the  continuum  flow  here;  the  amount  of  nonequilibrium  present  in  the  flow  seems  to 
be  adequately  handled  by  the  Navier-Stokes  equations  with  no-slip  boundary  condi¬ 
tions.  The  Mach  25  case  shows  an  increase  in  the  maximum  temperature  as  well  as 
a  slight  increase  in  the  shock  thickness;  however,  for  the  presently  considered  simple 
gas,  there  are  no  other  effects. 

Although  there  is  an  increasing  amount  of  nonequilibrium  in  the  flow  as  the 


(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  4.8:  Density  ratio  field  for  a  Mach  10  flow  of  argon  about  a  cylinder.  Density 
is  normalized  by  the  freestream  density.  The  normal  shock  relations 
predict  a  density  ratio  of  3.87. 


(a)  Kn  =  0.002  (b)  Kn  =  0.01 


(c)  Kn  =  0.05  (d)  Kn  =  0.25 


Figure  4.9:  Density  ratio  field  for  a  Mach  25  flow  of  argon  about  a  cylinder.  Density 
is  normalized  by  the  freestream  density.  The  normal  shock  relations 
predict  a  density  ratio  of  3.97. 
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(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  4. 10:  Temperature  field  for  a  Mach  10  flow  of  argon  about  a  cylinder.  The 
normal  shock  relations  predict  a  post-shock  temperature  of  6,460  K. 
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(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  4.11:  Temperature  field  for  a  Mach  25  flow  of  argon  about  a  cylinder.  The 
normal  shock  relations  predict  a  post-shock  temperature  of  39,400  K. 
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density  is  decreased,  the  temperature  fields  in  the  Kn  =  0.01  case  still  show  fairly 
good  agreement.  There  are,  however,  differences  in  the  shock  thickness,  with  DSMC 
predicting  a  more  diffuse  shock.  (The  Navier-Stokes  equations  are  known  to  poorly 
predict  the  interior  shock  structure  [25].)  Disagreement  is  also  seen  in  the  wake, 
where  the  continuum  hypothesis  is  expected  to  break  down  first.  The  shock  standoff 
distance  predicted  by  both  methods  is  nearly  the  same  as  is  the  thermal  boundary 
layer  thickness  at  the  stagnation  point. 

With  the  increase  in  nonequilibrium  present  in  the  Kn  =  0.05  case,  the  differences 
between  CFD  and  DSMC  temperature  field  predictions  are  more  pronounced,  with 
even  more  differences  for  the  higher  Mach  number.  Most  notable,  the  DSMC  shock 
is  much  thicker,  with  a  larger  stand-off  distance.  The  maximum  temperature  behind 
the  shock  has  also  decreased  from  the  higher  density  cases,  although  both  methods 
agree  fairly  well  in  their  predictions  of  the  peak  temperature  value  and  the  thermal 
boundary  layer  thickness  at  the  stagnation  point.  In  the  wake,  some  differences  are 
expected  with  DSMC  predicting  lower  temperatures. 

There  are  major  differences  present  for  the  Kn  =  0.25  case,  as  would  be  expected 
considering  the  low  density  of  the  flow.  There  is  a  further  decrease  in  the  magnitude 
of  the  peak  temperatures  behind  the  shock,  and  DSMC  predicts  a  higher  temperature 
than  CFD  for  the  Mach  25  case.  The  shock  thicknesses  and  stand-off  distances 
are  significantly  different,  as  is  the  thermal  boundary  layer  thickness,  with  DSMC 
predicting  a  more  diffusive  shock  and  boundary  layer. 

4.2.3  Stagnation  Line 

The  temperatures  predicted  along  the  stagnation  line  are  compared  in  Figures 
4.12  and  4.13,  again  for  Mach  10  and  Mach  25  respectively.  The  maximum  value  of 
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Khgll  is  also  plotted  against  the  right  axis. 

For  Kn  =  0.002,  all  of  the  CFD  simulations  agree  with  the  DSMC  results  (al¬ 
though  there  is  slight  disagreement  in  the  shock  for  the  Mach  25  case).  Figure  4.12(a) 
shows  a  typical  stagnation  line  temperature  distribution.  The  temperature  is  initially 
constant  at  the  free  stream  value  before  increasing  as  the  flow  is  compressed  in  the 
shock.  The  temperature  then  remains  fairly  constant  through  the  shock  layer  be¬ 
fore  decreasing  in  the  boundary  layer  to  the  wall  temperature.  Here,  the  only  areas 
where  there  is  a  significant  amount  of  continuum  breakdown  is  in  the  shock,  with  a 
smaller  amount  in  the  boundary  layer.  For  the  Kn  =  0.01  case,  there  is  still  good 
agreement  between  all  CFD  predictions  and  DSMC  for  the  post-shock  temperature 
distribution,  although  there  is  some  disagreement  in  the  shock  itself.  The  peak  value 
of  the  breakdown  parameter  is  also  higher  than  for  the  higher  density  case.  The 
Kn  =  0.05  case  shows  that  the  shock  layer  has  now  merged  with  the  boundary  layer. 
The  post-shock  temperatures  are  still  in  agreement,  although  there  are  differences 
among  the  CFD  cases.  The  Kn  =  0.25  case  shows  that  the  shock  and  the  boundary 
layer  have  merged  completely.  The  peak  temperature  is  lower,  and  there  are  signifi¬ 
cant  amounts  of  continuum  breakdown  from  the  shock  all  the  way  to  the  wall.  There 
is  also  some  significant  disagreement  between  each  of  the  CFD  cases. 

4.2.4  Surface  Properties 

The  surface  property  distributions  (pressure,  shear  stress  and  heat  flux)  for  each  of 
the  cases  is  examined  here.  The  surface  pressure,  in  the  form  of  a  pressure  coefficient, 
is  shown  in  Figures  4.14  and  4.15.  The  surface  pressure  is  the  least  sensitive  to 
nonequilibrium  of  all  the  surface  properties;  all  of  the  CFD  solutions  agree  well 
with  DSMC  for  all  but  the  most  rarefied  conditions  (Kn  =  0.25),  where  the  DSMC 
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Figure  4.12:  Temperature  profiles  along  the  stagnation  line  for  a  Mach  10  flow  of 
argon  about  a  cylinder.  The  maximum  value  of  KiiGLL,max  is  plotted  on 
the  right  axis.  Flow  is  from  left  to  right,  and  distance  is  normalized  by 
the  cylinder  radius. 
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Figure  4.13:  Temperature  profiles  along  the  stagnation  line  for  a  Mach  25  flow  of 
argon  about  a  cylinder.  The  maximum  value  of  Kiigll  is  plotted  on 
the  right  axis.  Flow  is  from  left  to  right;  distance  is  normalized  by  the 
cylinder  radius. 
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pressure  is  less  than  the  CFD  pressure  near  the  fore-body  (most  likely  due  to  the 
thinner  shock  predicted  by  CFD,  which  allows  the  CFD  solution  to  approach  the 
ideal  jump  condition  more  closely  than  the  DSMC  solution).  The  pressure  tensor  at 
this  Knudsen  number  is  also  most  likely  non-isotropic,  which  would  introduce  errors 
into  the  CFD  solution. 

Also  note  that  the  maximum  Kiigll  near  the  surface  is  plotted  against  the  right 
axis.  The  amount  of  continuum  breakdown  near  the  surface  is  highest  in  the  wake, 
as  expected,  and  the  amount  of  continuum  breakdown  increases  as  the  flow  becomes 
more  rarefied.  Nevertheless,  it  is  not  until  the  complete  merger  of  the  shock  and 
boundary  layer  in  the  Kn  =  0.25  case  that  there  is  an  effect  on  the  pressure  distribu¬ 
tion.  The  Kiigll  here  is  calculated  from  the  Gokgen  CFD  solution;  the  no-slip  CFD 
solution  shows  a  much  larger  amount  of  continuum  breakdown  near  the  surface  due 
to  the  larger  gradients  required  to  satisfy  the  no-slip  condition. 

The  shear  stress  is  the  most  sensitive  of  the  surface  properties  to  the  amount  of 
nonequilibrium  in  the  flow,  as  seen  in  Figures  4.16  and  4.17.  Notable  shear  stress 
differences  are  seen  in  the  Kn  =  0.01  case,  and  these  become  more  pronounced  as  the 
density  decreases.  The  best  agreement  with  the  DSMC  solution  is  achieved  with  the 
Gokgen  boundary  conditions;  this  most  certainly  is  the  reason  for  the  good  agreement 
between  the  Gokgen  CFD  and  DSMC  predictions  for  total  drag. 

The  heat  transfer  rate  distributions  are  shown  in  Figures  4.18  and  4.19.  Again,  the 
slip  boundary  conditions  improve  the  prediction  capabilities  of  CFD,  with  the  Gokgen 
condition  showing  the  most  agreement.  Again,  both  computational  methods  agree 
well  for  the  highest  density  cases  and  the  agreement  grows  worse  as  the  flow  becomes 
more  rarefied.  The  heat  transfer  rate  does  not  seem  to  be  sensitive  to  the  amount  of 


continuum  breakdown  near  the  surface  as  the  amount  of  error  between  DSMC  and 
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Figure  4.14:  Surface  pressure  coefficient  for  a  Mach  10  flow  of  argon  about  a  cylinder. 

The  maximum  value  of  KnQLL  near  the  surface  is  plotted  on  the  right 
axis.  $  is  the  angle  from  the  stagnation  point. 
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Figure  4.15:  Surface  pressure  coefficient  for  a  Mach  25  flow  of  argon  about  a  cylinder. 

The  maximum  value  of  KnQLL  near  the  surface  is  plotted  on  the  right 
axis.  $  is  the  angle  from  the  stagnation  point. 
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Figure  4.16:  Surface  friction  coefficient  for  a  Mach  10  flow  of  argon  about  a  cylinder. 

The  maximum  value  of  KnQLL  near  the  surface  is  plotted  on  the  right 
axis.  $  is  the  angle  from  the  stagnation  point. 
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Figure  4.17:  Surface  friction  coefficient  for  a  Mach  25  flow  of  argon  about  a  cylinder. 

The  maximum  value  of  KnQLL  near  the  surface  is  plotted  on  the  right 
axis.  $  is  the  angle  from  the  stagnation  point. 
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CFD  is  very  nearly  constant  along  the  entire  surface;  that  is,  the  disagreement  does 
not  necessarily  correlate  with  the  higher  values  of  Kiiqll  near  the  surface.  This  would 
suggest  that  errors  in  heat  transfer  rate  in  the  CFD  solutions  are  most  likely  caused 
by  differences  in  the  shock  thickness,  post-shock  peak  temperatures  and  thermal 
boundary  layer  thicknesses. 

4.2.5  Flow  Properties  Along  a  Line  at  <F  =  90° 

The  use  of  slip  boundary  conditions  has  been  shown  to  improve  the  agreement  in 
surface  property  prediction  between  DSMC  and  CFD.  In  this  section,  the  agreement 
in  flow  properties  is  investigated.  The  surface  properties  that  are  most  sensitive 
to  nonequilibrium  effects  are  the  shear  stress  and  heat  flux.  Therefore,  the  flow 
properties  that  are  examined  here  are  the  temperature  and  velocity. 

The  temperature  and  the  breakdown  parameter  values  are  plotted  along  a  line 
normal  to  the  surface  at  an  angle  of  90°  in  Figures  4.20  and  4.21,  where  the  distance 
from  the  wall  is  normalized  by  the  cylinder  radius.  A  trend  noted  here  is  that  the 
regions  where  values  of  Kiigll  are  highest  are  closest  to  the  surface  in  the  Knudsen 
layer,  which  is  where  a  temperature  jump  is  seen. 

Very  similar  trends  are  seen  with  the  velocity  magnitude,  which  is  plotted  along 
the  same  line  in  Figures  4.22  and  4.23.  Again,  there  are  significant  amount  of  non- 
equilibrium  near  the  wall,  with  a  velocity  jump  at  the  wall  for  the  highest  Knudsen 
number  cases.  Although  the  Gokgen  CFD  cases  show  the  best  agreement  in  sur¬ 
face  property  prediction,  they  do  not  necessarily  show  the  best  agreement  in  flow 
properties,  particularly  for  the  Kn  =  0.05  cases. 

For  the  Kn  =  0.002  cases,  there  is  very  good  agreement  between  all  the  CFD 
solutions  and  the  DSMC  solution,  although  there  is  a  modest  temperature  jump  near 
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Figure  4.18:  Surface  heating  coefficient  for  a  Mach  10  flow  of  argon  about  a  cylinder. 

The  maximum  value  of  KnQLL  near  the  surface  is  plotted  on  the  right 
axis.  is  the  angle  from  the  stagnation  point. 
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Figure  4.19:  Surface  heating  coefficient  for  a  Mach  25  flow  of  argon  about  a  cylinder. 

The  maximum  value  of  KnQLL  near  the  surface  is  plotted  on  the  right 
axis.  is  the  angle  from  the  stagnation  point. 
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Figure  4.20:  Temperature  along  a  line  normal  to  the  body  surface  at  $  =  90°  for  a 
Mach  10  flow  of  argon  about  a  cylinder.  Distance  is  normalized  by  the 
cylinder  radius  and  $  is  the  angle  from  the  stagnation  point.  Note  that 
the  wall  temperature  (at  S/R  =  0)  is  500  K. 
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Figure  4.21:  Temperature  along  a  line  normal  to  the  body  surface  at  $  =  90°  for  a 
Mach  25  flow  of  argon  about  a  cylinder.  Distance  is  normalized  by  the 
cylinder  radius  and  $  is  the  angle  from  the  stagnation  point.  Note  that 
the  wall  temperature  (at  S/R  =  0)  is  1500  K. 
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Figure  4.22:  Velocity  magnitude  along  a  line  normal  to  the  body  surface  at  4>  =  90° 
for  a  Mach  10  flow  of  argon  about  a  cylinder.  Distance  is  normalized  by 
the  cylinder  radius  and  $  is  the  angle  from  the  stagnation  point.  Note 
the  non-zero  velocity  slip  at  the  wall  (at  S/R  =  0). 
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Figure  4.23:  Velocity  magnitude  along  a  line  normal  to  the  body  surface  at  4>  =  90° 
for  a  Mach  25  flow  of  argon  about  a  cylinder.  Distance  is  normalized  by 
the  cylinder  radius  and  $  is  the  angle  from  the  stagnation  point.  Note 
the  non-zero  velocity  slip  at  the  wall  (at  S/R  =  0). 
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the  surface.  There  is  very  little  nonequilibrium  except  near  the  wall.  As  the  density 
decreases,  there  is  more  disagreement  among  the  several  CFD  solutions,  especially 
near  the  wall  but  the  agreement  improves  extending  out  into  the  flow  held.  For  the 
Kn  =  0.05,  Mach  10  case,  the  Maxwell  CFD  solution  agrees  best  with  DSMC,  while 
for  the  Mach  25  case,  the  Lockerby  solution  agrees  better  near  the  wall.  For  the 
Kn  =  0.25  case,  the  Gokgen  solution  shows  the  best  agreement. 

4.2.6  Slip  Quantities 

The  Gokgen  CFD  predictions  of  the  surface  properties  of  pressure,  shear  stress 
and  heat  flux  have  been  shown  to  agree  best  with  DSMC.  Here  the  velocity  slip  and 
temperature  jump  are  examined. 

The  velocity  slip  for  each  simulation  is  seen  in  Figures  4.24  and  4.25.  Although 
there  are  some  differences  in  the  actual  peak  velocity  slip  values  (especially  in  the 
Mach  25  cases),  the  qualitative  agreement  is  very  good.  Note  that  the  Lockerby 
boundary  conditions  agree  best  with  the  DSMC  data  because  of  its  use  of  the  cor¬ 
rect  velocity  slip  at  the  wall  rather  than  a  fictitious  slip  used  in  the  Maxwell  slip 
conditions,  as  discussed  in  Chapter  III.  The  Gokgen  solution  agrees  least;  it  was 
derived  specifically  to  match  the  wall  properties  of  shear  stress  and  heat  flux  at  the 
wall  rather  than  to  accurately  predict  the  velocity  slip  and  temperature  jump  at  the 
wall. 

The  temperature  jump  for  each  simulation  is  plotted  in  Figures  4.26  and  4.27. 
Here  there  is  less  agreement  between  CFD  and  DSMC,  although  the  qualitative 
agreement  at  Kn  =  0.002  is  fair.  Again,  the  Lockerby  boundary  conditions  show 
the  best  agreement;  DSMC  and  the  Gokgen  cases  differ  by  more  than  100%,  and  as 
much  as  500%,  at  the  stagnation  point  for  the  more  rarefied  cases. 
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Figure  4.24:  Velocity  slip  for  a  Mach  10  flow  of  argon  about  a  cylinder.  $  is  the 
angle  from  the  stagnation  point. 
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Figure  4.25:  Velocity  slip  for  a  Mach  25  flow  of  argon  about  a  cylinder.  $  is 
angle  from  the  stagnation  point. 
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Figure  4.26:  Temperature  jump  for  a  Mach  10  flow  of  argon  about  a  cylinder.  4>  is 
the  angle  from  the  stagnation  point. 


Temperature  Jump  [K]  Temperature  Jump  [K] 


-  DSMC 

- - CFD  (no-slip) 

- CFD  (Maxwell) 


(a)  Kn  =  0.002 


-  DSMC 

- - CFD  (no-slip) 


-  CFD  (Maxwell) 


-  DSMC 

- - CFD  (no-slip) 


-  CFD  (Maxwell) 


(b)  Kn  =  0.01 


-  DSMC 

- - CFD  (no-slip) 


(c)  Kn  =  0.05  (d)  Kn  =  0.25 

Figure  4.27:  Temperature  jump  for  a  Mach  25  flow  of  argon  about  a  cylinder.  $  is 
the  angle  from  the  stagnation  point. 
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4.2.7  Comparison  of  Solutions  Across  the  Knudsen  Layer 

Recall  that  the  disagreements  in  flow  properties  are  greatest  near  the  wall,  where 
there  is  more  nonequilibrium  than  elsewhere.  Here  the  simulations  are  analyzed  to 
determine  if  those  differences  are  actually  concentrated  in  the  Knudsen  layer.  The 
common  definition  of  Knudsen  layer  is  the  portion  of  the  flow  on  the  order  of  one 
mean  free-path  (MFP)  from  the  wall. 

The  analysis  looks  at  data  extracted  from  both  DSMC  and  CFD  solutions  in 
four  locations;  along  the  stagnation  line,  along  a  line  extending  normal  to  the  body 
at  a  point  45°  from  the  stagnation  point,  along  a  line  at  90°  from  the  stagnation 
point,  and  along  a  line  at  135°  degrees  (essentially  in  the  wake).  The  following  flow 
variables  are  examined:  velocity  magnitude,  temperature,  density  and  pressure.  The 
number  of  MFPs  from  the  wall  is  estimated  by  dividing  the  distance  from  the  wall 
by  the  local  MFP  (from  the  DSMC  solution).  This  is  only  an  estimate  clue  to  the 
variation  of  the  local  MFP  with  location.  The  percent  difference  is  obtained  by 
using  the  DSMC  value  as  the  base;  e.g.  the  percent  difference  in  temperature  is 
100  x  ( TCfd  —  Tdsmc) /Tdsmc- 

Figure  4.28  shows  the  results  for  temperature  difference  for  the  Mach  10,  Kn  =  0.01 
case.  Although  there  are  still  some  differences  between  CFD  and  DSMC  further  from 
the  wall,  the  locations  where  those  differences  are  greater  than  5%  is  concentrated 
in  the  region  within  10  MFPs  of  the  wall.  This  is  seen  most  clearly  in  the  Kn  =  0.01 
and  Kn  =  0.05  cases.  Note  that  the  Kn  =  0.25  flow  solutions  diverge  about  3  MFPs 
from  the  wall;  this  is  due  to  flow  differences  in  the  shock. 

Figure  4.29  shows  similar  results  for  the  Mach  25  cases.  Again,  the  regions  where 
the  differences  are  greater  than  5%  are  within  about  10  MFPs  of  the  wall.  For  this 
higher  velocity  case,  however,  the  Kn  =  0.25  solutions  do  not  reach  better  than  5% 
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Figure  4.28:  Percent  difference  in  temperature  predicted  in  Knudsen  layer  along  a 
line  normal  to  the  body  surface  at  <f>  =  90°  for  a  Mach  10  flow  of  argon 
about  a  cylinder.  Distance  is  normalized  by  the  cylinder  radius  and  <f> 
is  the  angle  from  the  stagnation  point.  For  Kn  =  0.25,  the  shock  starts 
about  3  MFPs  from  the  wall. 
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agreement  before  reaching  the  shock. 

Defining,  then,  the  Knudsen  layer  as  being  of  a  thickness  on  the  order  of  one 
MFP,  i.e.  it  is  10  MFPs  or  less,  the  differences  between  CFD  and  DSMC  are  indeed 
concentrated  in  the  Knudsen  layer.  This  statement,  however,  does  not  apply  near  the 
shock,  where  there  are  also  significant  differences  as  the  flow  becomes  more  rarefied. 

4.2.8  Free- Molecular  Flow  (Kn  — >  oo ) 

As  the  Knudsen  number  increases,  the  computational  results  for  surface  pressure, 
shear  stress  and  heat  transfer  should  approach  the  analytical  results  in  the  free- 
molecular,  or  collisionless,  limit  [31].  Figure  4.30  compares  the  non-dimensional 
surface  properties  obtained  from  the  Mach  10  flow  of  argon  DSMC  simulations  with 
those  obtained  analytically  for  free-molecular  flow.  As  a  further  comparison  of  the 
DSMC  implementation,  a  case  is  run  for  Kn  =  100,  which  is  essentially  collisionless. 
It  is  evident  that  as  the  Knudsen  number  increases,  the  DSMC  results  do  indeed 
approach  the  analytic  results.  In  fact,  for  Kn  =  100  the  DSMC  results  are  nearly 
identical  to  the  analytical  results. 

4.2.9  Computational  details 

Some  computational  details  of  the  simulations  discussed  here  are  given  in  Table 
4.5.  It  should  be  noted  that  the  Gokgen  CFD  cases  generally  took  longer  to  converge 


than  the  other  CFD  cases. 
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Figure  4.29:  Percent  difference  in  temperature  predicted  in  Knudsen  layer  along  a 
line  normal  to  the  body  surface  at  <f>  =  90°  for  a  Mach  25  flow  of  argon 
about  a  cylinder.  Distance  is  normalized  by  the  cylinder  radius  and  <f> 
is  the  angle  from  the  stagnation  point.  For  Kn  =  0.25,  the  shock  starts 
about  1.5  MFPs  from  the  wall. 
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Figure  4.30:  Comparisons  of  surface  properties  predicted  by  DSMC  from  the  contin¬ 
uum  to  the  free-molecular  regimes.  As  the  Knudsen  number  increases, 
the  DSMC  solutions  tend  toward  the  free-molecular  analytical  limits. 
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Tabic  4.5:  Computational  details  for  a  flow  of  argon  about  a  cylinder.  Total  CPU 
time  is  the  wall  time  multiplied  by  the  number  of  CPUs. 

Mach  10 


Knoo 

DSMC 

Cells 

Particles 

Time  Steps 

CPUs 

Total  CPU  Time* [hours 

0.002 

394,250 

287.0x10® 

1,050,000 

128 

57,186 

0.01 

68,497 

26.8x10® 

430,000 

32 

1,827 

0.05 

18,818 

7.1x10® 

430,000 

32 

499 

0.25 

24,452 

8.2x10® 

290,000 

32 

387 

Knoo 

CFD 

Cells 

Iterations 

CPUs 

Total  CPU  Time* [hours 

0.002 

80,000 

25,000 

8 

150 

0.01 

40,000 

25,000 

8 

75 

0.05 

40,000 

25,000 

8 

75 

0.25 

40,000 

25,000 

8 

75 

Mach  25 


Knoo 

DSMC 

Cells 

Particles 

Time  Steps 

CPUs 

Total  CPU  Time* [hours 

0.002 

102,626 

75.1x10® 

420,000 

64 

5,047 

0.01 

71,997 

37.7x10® 

300,000 

32 

1,793 

0.05 

22,523 

13.1x10® 

200,000 

16 

429 

0.25 

18,395 

6.2x10® 

200,000 

8 

198 

Knoo 

CFD 

Cells 

Iterations 

CPUs 

Total  CPU  Time* [hours 

0.002 

80,000 

25,000 

8 

150 

0.01 

40,000 

25,000 

8 

75 

0.05 

40,000 

25,000 

8 

75 

0.25 

40,000 

25,000 

8 

75 

Approximate 
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4.3  Nitrogen 

The  case  of  a  hypersonic  flow  of  nitrogen  about  a  cylinder  is  now  considered. 
Again,  DSMC  and  CFD  simulations  are  compared. 

All  temperature  results  from  the  DSMC  simulations  are  obtained  by  averaging 
the  translational  and  rotational  temperatures,  as  shown  in  Eq.  3.20  unless  otherwise 
noted. 

Tables  4.6  and  4.7  compare  the  total  drag  and  the  peak  heat  transfer  rate  pre¬ 
dicted  by  both  computational  methods.  The  differences  between  DSMC  and  CFD 
are  also  shown  graphically  in  Figures  4.32  and  4.33.  Here,  the  peak  heat  transfer  rate 
is  obtained  by  averaging  over  the  surface  area  within  one  degree  of  the  stagnation 
point.  For  CFD,  these  quantities  are  calculated  for  the  no-slip  boundary  conditions; 
the  Gokgen  slip  conditions  [29],  which  gave  the  best  agreement  with  DSMC  for  the 
surface  properties;  and  the  Lockerby  wall-function  slip  conditions  [49],  which  gave 
the  best  agreement  with  DSMC  for  the  slip  data. 

Figure  4.31  illustrates  the  percentage  of  total  drag  due  to  pressure  and  friction 
forces,  for  both  DSMC  and  CFD.  As  with  argon,  the  percentage  of  total  drag  due  to 
friction  increases  from  about  5%  at  Kn  =  0.002  to  about  20%  at  Kn  =  0.25  and  the 
vast  majority  of  the  drag  is  due  to  pressure  forces.  That  the  difference  in  predicted 
total  drag  between  CFD  and  DSMC  is  also  due  mostly  to  the  differences  predicted 
in  the  friction  forces,  as  shown  in  Figure  4.32. 

Although  here  a  gas  with  internal  degrees  of  freedom  is  considered,  the  additional 
possibility  of  thermodynamic  nonequilibrium  does  not  seem  to  significantly  affect  the 
surface  properties,  or  the  total  drag  and  stagnation  heat  transfer.  As  with  the  flow  of 
argon  noted  previously,  the  slip  boundary  conditions  improve  the  agreement  between 
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Table  4.6:  Total  drag  for  a  flow  of  nitrogen  about  a  cylinder. 

Mach  10 


Drag/Length  [N/m]  (%  Difference) 

Klloo 

DSMC 

CFD 

no-slip 

Gokgen 

Lockerby 

0.002 

162.4 

162.4  (0.0%) 

162.2  (-0.2%) 

162.4  (0.0%) 

0.01 

34.23 

34.37  (0.4%) 

34.17  (-0.2%) 

34.32  (0.3%) 

0.05 

7.346 

7.780  (5.9%) 

7.456  (1.5%) 

7.603  (3.5%) 

0.25 

1.714 

2.056  (19.9%) 

1.702  (-0.7%) 

1.868  (9.0%) 

Mach  25 


Drag/Length  [N/m]  (%  Difference) 

Klloo 

DSMC 

CFD 

no-slip 

Gokgen 

Lockerby 

0.002 

1020. 

1022.  (0.2%) 

1021.  (0.1%) 

1022.  (0.2%) 

0.01 
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Figure  4.31:  Percentage  of  total  drag  due  to  pressure  and  friction  forces  for  a  flow  of 
nitrogen  about  a  cylinder. 
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Table  4.7:  Peak  heat  transfer  rate  for  a  flow  of  nitrogen  about  a  cylinder. 

Mach  10 
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Figure  4.32:  Total  drag  difference  from  DSMC  for  flow  of  nitrogen  about  a  cylinder. 

Note  that  there  is  a  negative  contribution  to  total  drag  due  to  friction 
for  the  Kn  =  0.25,  Mach  25  case,  bringing  the  total  drag  difference  down 
to  2.8%. 
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Figure  4.33:  Peak  heat  transfer  rate  difference  from  DSMC  for  flow  of  nitrogen  about 
a  cylinder. 

the  two  simulation  methods,  with  the  Gokgen  slip  boundary  conditions  showing  the 
best  agreement.  In  some  cases,  there  is  better  agreement  with  the  nitrogen  cases 
than  with  the  argon  cases;  for  example,  for  the  Mach  10,  Kn  =  0.25  case,  the  total 
drag  computed  using  the  Gokgen  boundary  conditions  is  within  3%  of  the  DSMC 
solution  with  nitrogen,  although  it  is  just  over  5%  when  using  argon. 

The  peak  heat  transfer  rates  also  show  similar  trends  for  nitrogen  as  with  argon, 
although  here  the  largest  difference  seen  in  the  nitrogen  results  using  the  Gokgen 
boundary  conditions  is  8%,  compared  to  about  6%  with  argon.  The  Mach  25  heat 
transfer  results  show  a  larger  difference  when  considering  nitrogen.  The  relatively 
large  differences  for  the  Mach  25,  Kn  =  0.002  case  are  most  likely  due  to  computa¬ 
tional  artifacts,  rather  than  complexity  due  to  the  nonequilibrium  nature  of  the  flow, 
as  will  be  discussed  below. 

Overall,  however,  the  current  results  for  nitrogen  are  very  similar  to  those  for 
argon.  The  Gokgen  conditions  give  the  best  results,  with  total  drag  being  within  3% 
of  the  DSMC  results,  and  the  peak  heating  rates  being  within  about  5%. 
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Again,  given  that  the  Gokgen  condition  shows  the  best  agreement  with  the  DSMC 
solutions,  the  CFD  field  results  shown  below  are  taken  from  the  Gdkgen  CFD  case. 

4.3.1  Continuum  Breakdown 

The  breakdown  parameter  is  again  calculated  using  both  the  CFD  and  the  DSMC 
solutions.  As  before,  the  only  causes  of  breakdown  to  the  continuum  hypothesis 
expected  here  are  in  regions  of  high  gradients  (such  as  the  shock  and  boundary 
layer)  and  regions  of  rarefaction  (such  as  in  the  wake).  Thermal  nonequilibrium  is 
expected  to  be  higher  in  regions  where  the  breakdown  parameter  values  are  highest; 
thus,  the  nonequilibrium  between  translational  and  rotational  modes  is  expected 
in  the  shock  wave,  and  the  wake.  As  the  flow  becomes  more  rarefied  in  the  higher 
Knudsen- number  cases,  the  degree  of  thermal  nonequilibrium  is  expected  to  increase. 

Figures  4.34  and  4.35  illustrate  the  amount  of  continuum  breakdown  in  the  flow  as 
it  becomes  more  rarefied.  As  with  the  argon  flow,  the  maximum  gradient  length  local 
Knudsen  number  is  computed  from  the  DSMC  (top)  and  CFD  (bottom)  solutions. 
Light  gray  regions  correspond  to  Kiigll  <  0.05,  dark  gray  regions  correspond  to 
Kugll  <  0.10  and  black  regions  correspond  to  KnGLL  >  0.10.  Continuum  breakdown 
is  expected  in  the  same  regions  of  the  flow,  namely  in  the  shock  region,  in  a  thin 
boundary  layer  along  the  surface  and  in  a  region  of  flow  expansion  around  the  top  of 
the  cylinder.  The  additional  thermal  nonequilibrium  present  due  to  internal  degrees 
of  freedom  does  not  affect  the  amount  of  continuum  breakdown  predicted  by  the 
breakdown  parameter. 

4.3.2  Flow  Field  Properties 

The  density  ratio  fields,  where  the  density  is  normalized  by  the  freestream  density, 
can  be  seen  in  Figures  4.36  and  4.37.  Note  that  the  maximum  density  ratio  behind  the 


(c)  Kn  =  0.05  (d)  Kn  =  0.25 

Figure  4.34:  Kiigll  field  for  a  Mach  10  flow  of  nitrogen  about  a  cylinder.  Light  gray 


regions  correspond  to  Kiigll  <  0.05,  dark  gray  regions  correspond  to 
Kiigll  <  0.10  and  black  regions  correspond  to  Kiigll  >  0.10. 
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(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  4.35:  Kiigll  field  for  a  Mach  25  flow  of  nitrogen  about  a  cylinder.  Light  gray 
regions  correspond  to  Kiigll  <  0.05,  dark  gray  regions  correspond  to 
Kiigll  <  0.10  and  black  regions  correspond  to  Kiigll  >  0.10. 
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DSMC 

(d)  Kn  =  0.25 


Figure  4.36:  Density  ratio  field  for  a  Mach  10  flow  of  nitrogen  about  a  cylinder.  The 
normal  shock  relations  predict  a  density  ratio  of  5.71. 


shock  is  about  6,  agreeing  with  the  theoretical  limit  for  a  diatomic  gas  in  equilibrium 
and  the  results  of  the  normal  shock  relations,  p/po  =  5.87  (Mach  10)  and  p/po  =  5.95 
(Mach  25).  Overall  agreement  between  CFD  and  DSMC  is  good,  with  some  small 
differences  in  the  shock  structure  at  the  higher  Knudsen  number  conditions.  As  with 
argon,  as  the  freestream  density  decreases,  the  shock  becomes  much  more  diffuse. 

The  translational/rotational  temperature  fields  are  seen  in  Figures  4.38  and  4.39. 
The  temperature  results  for  the  Mach  10  flows  are  very  similar  to  those  seen  above 
with  the  argon  flow.  However,  the  Mach  25  results  exhibit  some  significant  differ- 
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ences.  Although  the  CFD  and  DSMC  methods  agree  well  for  Kn  =  0.002,  the  peak 
temperature  seen  in  the  narrow  band  immediately  behind  the  shock  in  Figure  4.39(a) 
shows  the  effects  of  vibrational  energy  excitation;  that  is,  as  the  translational  and 
rotational  temperatures  rise,  energy  is  transferred  into  the  vibrational  mode,  decreas¬ 
ing  the  translational  temperature.  This  will  be  seen  more  clearly  in  the  temperature 
profiles  along  the  stagnation  line  that  are  examined  below.  The  Kn  =  0.01  cases 
also  show  some  effect  of  the  vibrational  energy  excitation,  and  the  CFD  solution  no 
longer  agrees  with  the  DSMC  solution.  The  Kn  =  0.05  and  Kn  =  0.25  cases  again 
show  large  differences,  with  the  DSMC  shocks  being  much  more  diffuse,  the  DSMC 
thermal  boundary  layer  being  much  thicker,  and  the  DSMC  temperature  in  the  wake 
being  lower.  In  each  of  these  cases  the  differences  are  more  extreme  than  was  noted 
for  the  argon  flow. 

The  vibrational  temperature  fields  can  be  seen  in  Figures  4.40  and  4.41.  Here 
it  can  be  seen  that,  for  the  Mach  10  cases,  the  amount  of  vibrational  excitation  is 
minimal.  The  peak  translational/rotational  temperatures  seen  in  the  Kn  =  0.002 
case  is  around  4000  K,  while  the  vibrational  temperature  peaks  at  about  1200  K. 
While  CFD  and  DSMC  are  in  good  agreement  for  the  vibrational  temperature  field 
in  the  fore-body  region,  CFD  predicts  higher  vibrational  temperatures  in  the  wake. 
As  the  flow  becomes  more  rarefied,  the  level  of  vibrational  excitation  decreases  as 
there  are  fewer  collisions  with  which  to  equilibrate  the  vibrational  and  translational 
energy  modes.  For  the  Kn  =  0.05  and  Kn  =  0.25  cases,  the  vibrational  temperature 
does  not  exceed  the  wall  temperature  value  of  500  K,  and  most,  if  not  all,  of  the 
vibrational  excitation  is  due  to  the  wall  boundary  condition. 

The  Mach  25  cases  show  a  larger  amount  of  vibrational  excitation,  as  expected, 
due  to  the  much  higher  temperatures.  The  peak  translational/rotational  temperature 


105 


(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  4.38:  Translational/rotational  temperature  field  for  a  Mach  10  flow  of  nitro¬ 
gen  about  a  cylinder.  The  normal  shock  relations  predict  a  post-shock 
temperature  of  4,078  K. 
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(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  4.39:  Translational/rotational  temperature  field  for  a  Mach  25  flow  of  nitro¬ 
gen  about  a  cylinder.  The  effects  of  vibrational  nonequilibrium  can 
be  seen  in  the  narrow  band  of  high  temperature  directly  behind  the 
gas,  as  energy  is  transferred  from  the  translational/rotational  modes 
to  the  vibrational  mode.  The  normal  shock  relations  give  a  post-shock 
temperature  of  24,500  K,  which  is  slightly  higher  than  the  maximum 
temperature  in  the  simulations. 
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(c)  Kn  =  0.05  (d)  Kn  =  0.25 

Figure  4.40:  Vibrational  temperature  field  for  a  Mach  10  flow  of  nitrogen  about 


a  cylinder.  The  vibrational  energy  modes  are  activated  in  significant 
amounts  only  for  the  lower  Knudsen  number  cases.  For  the  higher 
Knudsen  number  cases,  the  vibrational  modes  are  activated  mostly  by 
the  cylinder  surface,  where  the  vibrational  temperature  is  500  K. 
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(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  4.41:  Vibrational  temperature  field  for  a  Mach  25  flow  of  nitrogen  about  a 
cylinder.  The  vibrational  energy  modes  are  much  more  significant  than 
at  Mach  10  due  to  the  higher  temperatures  at  Mach  25.  Note  that  the 
vibrational  temperature  at  the  wall  is  1500  K. 
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for  the  Kn  =  0.002  case  is  around  20,000  K,  and  the  peak  vibrational  temperature  is 
very  nearly  the  same;  thus,  there  are  sufficient  collisions  such  that  the  different  energy 
modes  very  nearly  reach  equilibrium.  Again,  the  two  computational  methods  agree 
well  in  their  predictions  for  vibrational  temperature  for  this  and  the  Kn  =  0.01  case. 
The  CFD  solutions  significantly  under-predict  the  amount  of  vibrational  excitation 
for  the  more  rarefied  cases.  DSMC  predicts  a  significant  amount  of  vibrational 
excitation  within  the  shock  for  the  Kn  =  0.25  case,  while  that  in  the  CFD  solution 
is  mostly  due  to  the  wall  boundary  conditions. 

4.3.3  Stagnation  Line  Properties 

The  temperature  profiles  along  the  stagnation  line  for  the  Mach  10  cases  are 
shown  in  Figure  4.42.  Here,  the  translational,  rotational,  vibrational  temperatures 
and  the  averaged  translational/rotational  temperature  from  the  DSMC  solutions  are 
compared  with  the  translational/rotational  and  vibrational  temperatures  from  the 
Gokgen  CFD  solutions.  The  Kn  =  0.002  case  shows  that  there  is  very  little  none¬ 
quilibrium  except  in  the  shock.  In  the  shock,  there  is  a  slight  peak  in  translational 
temperature  prior  to  equilibrating  with  the  rotational  mode.  (This  effect  is  seen  more 
clearly  in  the  Kn  =  0.01  case.)  The  averaged  DSMC  temperature  and  the  CFD  tem¬ 
perature  agree  very  well,  although  there  is  a  small  amount  of  translational-rotational 
nonequilibrium.  There  is  also  a  significant  amount  of  vibrational  nonequilibrium  as 
the  temperatures  are  not  high  enough  to  significantly  excite  the  vibrational  mode. 
The  two-temperature  model  of  LeMANS  is  able  to  model  this  nonequilibrium  very 
satisfactorily. 

As  the  flow  becomes  more  rarefied,  the  amount  of  translational/rotational  none¬ 
quilibrium  in  the  shock  increases.  For  the  Kn  =  0.01  case,  however,  these  two  energy 
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(c)  Kn  =  0.05 
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Figure  4.42:  Temperature  profiles  along  the  stagnation  line  for  a  Mach  10  flow  of 
nitrogen  about  a  cylinder.  Flow  is  from  left  to  right,  with  the  distance 
normalized  by  the  cylinder  radius.  The  translational  and  rotational 
temperatures  from  DSMC  are  shown  along  with  the  averaged  transla¬ 
tional  / rotational  temperature. 
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(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  4.43:  Temperature  profiles  along  the  stagnation  line  for  a  Mach  25  flow  of 
nitrogen  about  a  cylinder.  Flow  is  from  left  to  right  with  the  distance 
normalized  by  the  cylinder  radius. 


Maximum  Knril  Maximum  Kn. 
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modes  are  quickly  brought  back  into  equilibrium  after  the  shock.  For  the  Kn  =  0.05 
and  Kn  =  0.25  cases,  however,  there  are  too  few  collisions  for  the  translational  and 
rotational  modes  to  equilibrate  prior  to  reaching  the  boundary  layer  (for  Kn  =  0.05) 
or  the  wall  (for  Kn  =  0.25).  The  vibrational  temperature  profiles  for  these  last  two 
cases  strongly  imply  that  most,  if  not  all,  of  the  vibrational  excitation  is  due  to  the 
wall  boundary  conditions. 

The  Mach  25  temperature  profiles  are  seen  in  Figure  4.42.  The  additional  vibra¬ 
tional  excitation  clue  to  the  high  temperatures  in  the  post-shock  region  is  clearly 
shown.  There  is  additional  translational-rotational  nonequilibrium,  even  in  the 
Kn  =  0.002  case.  However,  all  temperatures,  including  the  vibrational  temperature, 
equilibrate  prior  to  the  boundary  layer.  Similar  trends  as  with  Mach  10  are  noted 
for  the  more  rarefied  cases.  Note,  however,  that  CFD  under-predicts  the  vibrational 
temperature  in  the  Kn  =  0.05  and  Kn  =  0.25  cases. 

4.3.4  Surface  Properties 

The  surface  properties  for  the  flow  of  nitrogen  about  the  cylinder  are  examined  in 
this  section.  Interestingly,  there  is  no  significant  difference  between  the  trends  seen 
here  and  those  noted  for  the  argon  cases,  notwithstanding  the  additional  thermal 
nonequilibrium  present  in  the  flows.  At  most,  the  surface  property  predictions  (in 
terms  of  the  total  drag  and  peak  heat  transfer  rate)  differ  by  a  few  percentage  points 
more  than  they  did  for  the  simple  gas  cases. 

The  surface  pressure  distributions,  shown  in  Figures  4.44  and  4.45  are  very  similar 
to  those  for  argon.  In  particular,  all  of  the  simulations  agree  very  well,  except  for 
the  Kn  =  0.25  cases,  where  CFD  predicts  a  higher  pressure  than  DSMC. 

The  surface  shear  stress  distributions  are  seen  in  Figures  4.46  and  4.47.  Again, 
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(a)  Kn  =  0.002 


(b)  Kn  =  0.01 
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(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  4.44:  Surface  pressure  coefficient  for  a  Mach  10  flow  of  nitrogen  about  a  cylin¬ 
der.  The  maximum  Kiiqll  near  the  surface  is  plotted  on  the  right  axis. 
<3>  is  the  angle  from  the  stagnation  point. 
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(a)  Kn  =  0.002 


(b)  Kn  =  0.01 
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Figure  4.45:  Surface  pressure  coefficient  for  a  Mach  25  flow  of  nitrogen  about  a  cylin¬ 
der.  The  maximum  Kiiqll  near  the  surface  is  plotted  on  the  right  axis. 
<3>  is  the  angle  from  the  stagnation  point. 
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the  results  are  very  similar  to  those  for  argon,  although  CFD  diverges  more  from 
DSMC  as  the  Knudsen  number  increases.  Nevertheless,  the  slip  boundary  conditions 
greatly  improve  the  agreement,  with  the  Gokgen  boundary  conditions  showing  the 
best  agreement  with  DSMC. 

The  heat  transfer  rate  distributions  are  shown  in  Figures  4.48  and  4.49.  Here,  the 
Kn  =  0.002  cases  the  current  differences  of  about  3%  are  reasonable.  These  results 
are  also  very  similar  the  argon  cases,  with  the  CFD  results  diverging  from  DSMC  at 
the  higher  Knudsen  number  conditions  and  with  the  Tye  2  slip  conditions  giving  the 
best  agreement. 

For  Kn  =  0.01,  the  total  drag  predicted  by  CFD  is  still  within  1%  of  that  predicted 
by  DSMC,  for  Mach  10,  and  within  2%  for  Mach  25,  as  shown  in  Table  4.6.  The 
peak  heating  also  differs  by  about  2%  for  all  cases. 

The  total  drag  for  the  Kn  =  0.05  case  solutions  are  within  4%  of  the  DSMC  case, 
again  with  the  Gokgen  boundary  condition  showing  the  best  agreement.  While  the 
no-slip  case  predicts  a  peak  heat  transfer  rate  about  9%  higher  than  DSMC,  the  slip 
cases  show  better  agreement,  with  the  Gokgen  case  being  within  5%. 

For  the  most  rarefied  case  considered,  the  peak  heat  transfer  rates  differ  by  as 
much  as  almost  30%  for  the  no-slip  case  at  Mach  25,  but  the  Gokgen  case  shows  the 
best  agreement  with  just  over  5%  difference  at  Mach  10  and  8%  difference  at  Mach 
25. 

4.3.5  Flow  Properties  Along  a  Line  at  <f>  =  90° 

The  solution  temperatures  (including  the  translational,  rotational,  and  trans¬ 
lational/rotational  averaged  temperature)  are  plotted  along  a  line  normal  to  the 
surface  at  an  angle  of  90°  as  shown  in  Figures  4.50  and  4.51.  Note  the  increasing 
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(b)  Kn  =  0.01 
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Figure  4.46:  Surface  friction  coefficient  for  a  Mach  10  flow  of  nitrogen  about  a  cylin¬ 
der.  The  maximum  Kiiqll  near  the  surface  is  plotted  on  the  right  axis. 
<3>  is  the  angle  from  the  stagnation  point. 
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Figure  4.47:  Surface  friction  coefficient  for  a  Mach  25  flow  of  nitrogen  about  a  cylin¬ 
der.  The  maximum  Kiiqll  near  the  surface  is  plotted  on  the  right  axis. 
4>  is  the  angle  from  the  stagnation  point. 
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Figure  4.48:  Surface  heating  coefficient  for  a  Mach  10  flow  of  nitrogen  about  a  cylin¬ 
der.  The  maximum  Kiigll  near  the  surface  is  plotted  on  the  right  axis. 
4>  is  the  angle  from  the  stagnation  point. 


Maximum  Knril  Maximum  Kn. 


Heating  Coefficient,  CH  Heating  Coefficient,  C 


119 


-  DSMC 

- CFD  (no-slip) 


* 

E 


-  DSMC 

- - CFD  (no-slip) 


(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


-  DSMC 

- CFD  (no-slip) 


-  DSMC 

- - CFD  (no-slip) 


(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  4.49:  Surface  heating  coefficient  for  a  Mach  25  flow  of  nitrogen  about  a  cylin¬ 
der.  The  maximum  Kiigll  near  the  surface  is  plotted  on  the  right  axis, 
is  the  angle  from  the  stagnation  point. 
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amount  of  nonequilibrium  as  quantified  by  values  of  Kiiqll  as  the  Knudsen  number 
increases.  Although  there  is  some  translational-rotational  nonequilibrium,  there  is 
reasonable  agreement  between  the  Gokgen  and  Lockerby  solutions  and  DSMC,  with 
the  Lockerby  solutions  being  slightly  better. 

There  is  also  good  agreement  between  the  CFD  solutions  and  the  DSMC  solution 
for  vibrational  temperature,  shown  in  Figures  4.52  and  4.53  for  most  of  the  cases. 
The  Mach  10,  Kn  =  0.25  case  shows  good  agreement  between  DSMC  and  the  Gokgen 
CFD  solution,  and  the  Mach  25,  Kn  =  0.25  case  shows  a  very  large  disagreement 
between  DSMC  and  all  of  the  CFD  solutions.  In  this  case,  DSMC  shows  a  significant 
amount  of  vibrational  excitation  away  from  the  wall,  while  CFD  predicts  vibrational 
excitation  only  clue  to  the  wall  boundary  conditions. 

Again,  very  similar  trends  are  seen  with  the  velocity  magnitude,  which  is  plotted 
along  the  same  line  in  Figures  4.54  and  4.55.  Here,  the  Gokgen  and  Lockerby  so¬ 
lutions  show  reasonable  agreement  for  most  of  the  cases,  with  the  Gokgen  solution 
being  better  at  Kn  =  0.25. 

4.3.6  Slip  Quantities 

The  velocity  slip  for  each  nitrogen  simulation  is  seen  in  Figures  4.56  and  4.57. 
Again,  the  results  are  very  similar  to  those  obtained  for  argon.  Although  there  are 
some  differences  in  the  actual  peak  velocity  slip  values,  the  qualitative  agreement  is 
very  good.  Again,  note  that  the  Lockerby  boundary  conditions  agree  best  with  the 
DSMC  data  because  of  its  use  of  the  correct  velocity  slip  at  the  wall  rather  than 
a  fictitious  slip,  as  discussed  in  Chapter  III.  Recall  that  the  Gokgen  solution  was 
derived  specifically  to  match  the  wall  properties  of  shear  stress  and  heat  flux  at  the 
wall  rather  than  to  accurately  predict  the  velocity  slip  and  temperature  jump  at  the 
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(c)  Kn  =  0.05  (d)  Kn  =  0.25 

Figure  4.50:  Translational/rotational  temperatures  along  a  line  normal  to  the  body 
surface  at  $  =  90°  for  a  Mach  10  flow  of  nitrogen  about  a  cylinder. 
Distance  is  normalized  by  the  cylinder  radius  and  $  is  the  angle  from 
the  stagnation  point.  The  separate  translational  and  rotational  temper¬ 
atures  from  the  DSMC  simulations  are  shown  along  with  the  weighted- 
average  temperature.  Note  that  the  wall  temperature  (at  S/R  =  0)  is 
500  K. 
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Figure  4.51:  Translational/rotational  temperatures  along  a  line  normal  to  the  body 
surface  at  $  =  90°  for  a  Mach  25  flow  of  nitrogen  about  a  cylinder. 
Distance  is  normalized  by  the  cylinder  radius  and  $  is  the  angle  from 
the  stagnation  point.  The  separate  translational  and  rotational  temper¬ 
atures  from  the  DSMC  simulations  are  shown  along  with  the  weighted 
average  temperature.  Note  that  the  wall  temperature  (at  S/R  —  0)  is 
1500  K. 
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Figure  4.52:  Vibrational  temperature  along  a  line  normal  to  the  body  surface  at 
4>  =  90°  for  a  Mach  10  flow  of  nitrogen  about  a  cylinder.  Distance  is 
normalized  by  the  cylinder  radius  and  4>  is  the  angle  from  the  stagnation 
point.  Vibrational  temperature  decreases  away  from  the  wall  for  the 
higher  Knudsen  number;  vibrational  activation  is  due  only  to  the  wall 
boundary  conditions,  where  the  vibrational  temperature  is  500  K. 
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Figure  4.53:  Vibrational  temperature  along  a  line  normal  to  the  body  surface  at 
<3>  =  90°  for  a  Mach  25  flow  of  nitrogen  about  a  cylinder.  Distance  is 
normalized  by  the  cylinder  radius  and  <f>  is  the  angle  from  the  stagnation 
point.  None  of  the  CFD  solutions  agree  well  with  DSMC  at  Kn  =  0.25; 
the  vibrational  temperature  away  from  the  wall  in  the  CFD  cases  are 
near  the  wall  temperature  of  1500  K,  while  DSMC  predicts  increased 
vibrational  activation  further  from  the  wall. 
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Figure  4.54:  Velocity  magnitude  along  a  line  normal  to  the  body  surface  at  <f>  =  90° 
for  a  Mach  10  flow  of  nitrogen  about  a  cylinder.  Distance  is  normalized 
by  the  cylinder  radius  and  4>  is  the  angle  from  the  stagnation  point. 
Note  the  non- zero  velocity  slip  at  the  wall  (at  S/R  —  0). 
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Figure  4.55:  Velocity  magnitude  along  a  line  normal  to  the  body  surface  at  4>  =  90° 
for  a  Mach  25  flow  of  nitrogen  about  a  cylinder.  Distance  is  normalized 
by  the  cylinder  radius  and  4>  is  the  angle  from  the  stagnation  point. 
Note  the  non- zero  velocity  slip  at  the  wall  (at  S/R  —  0). 
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wall. 

The  temperature  jump  is  plotted  in  Figures  4.58  and  4.59.  Here,  the  DSMC 
translational/rotational  averaged  temperature  jump  is  shown  along  with  the  actual 
translational  and  rotational  temperature  jumps  (with  the  exception  of  the  Kn  =  0.002 
case,  for  which  the  rotational  temperature  jump  data  is  unavailable),  and  compared 
to  the  CFD  temperature  jump.  It  is  interesting  to  note  the  amount  of  translational- 
rotational  nonequilibrium  present  at  the  surface  as  seen  in  the  differences  between 
the  translational  and  rotational  temperature  jumps.  Again,  there  is  fair  agreement 
at  the  lower  Knudsen  numbers.  The  differences  are  very  pronounced  at  the  highest 
Knudsen  numbers,  with  the  Gokgen  cases  showing  a  nearly  500%  difference  at  the 
stagnation  point,  and  the  Lockerby  cases  show  the  best  agreement. 

The  vibrational  temperature  jump  for  each  simulation  (with  the  exception  of  the 
Kn  =  0.002  case,  for  which  the  vibrational  temperature  jump  is  not  available)  is 
plotted  in  Figures  4.60  and  4.61.  Qualitatively,  the  vibrational  temperature  jump 
curves  look  similar  to  the  translational/rotational  temperature  jump  curves. 

For  the  Mach  10,  Kn  =  0.01  and  Kn  =  0.05  cases,  there  is  reasonable  agree¬ 
ment  between  CFD  and  DSMC,  notwithstanding  the  statistical  scatter  present  in 
the  DSMC  data.  The  large  scatter  is  due  to  the  small  number  of  particles  with  a 
significant  amount  of  vibrational  energy,  which  is  expected  in  a  flow  with  the  rela¬ 
tively  low  temperatures.  For  the  Mach  10,  Kn  =  0.25  case,  the  Lockerby  CFD  case 
agrees  well  in  the  fore-body  region,  but  predicts  a  positive  vibrational  temperature 
jump  in  the  wake  while  DSMC  predicts  a  negative  vibrational  temperature  jump. 
Incidently,  the  negative  vibrational  temperature  jump  also  supports  the  conclusion 
that  the  vibrational  energy  present  in  the  flow  at  the  most  rarefied  conditions  is  due 
to  the  wall  boundary  conditions. 
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Figure  4.56:  Velocity  slip  for  a  Mach  10  flow  of  nitrogen  about  a  cylinder.  $  is  the 
angle  from  the  stagnation  point. 
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Figure  4.57:  Velocity  slip  for  a  Mach  25  flow  of  nitrogen  about  a  cylinder.  $  is  the 
angle  from  the  stagnation  point. 
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Figure  4.58:  Translational/rotational  temperature  jump  for  a  Mach  10  flow  of  nitro¬ 
gen  about  a  cylinder.  Separate  translational  and  rotational  temperature 
jumps  from  DSMC  are  plotted  as  well  as  the  average  temperature  jump. 
4>  is  the  angle  from  the  stagnation  point. 
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Figure  4.59:  Translational/rotational  temperature  jump  for  a  Mach  25  flow  of  nitro¬ 
gen  about  a  cylinder.  Separate  translational  and  rotational  temperature 
jumps  from  DSMC  are  plotted  as  well  as  the  average  temperature  jump. 
4>  is  the  angle  from  the  stagnation  point. 
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Figure  4.60:  Vibrational  temperature  jump  for  a  Mach  10  flow  of  nitrogen  about  a 
cylinder.  $  is  the  angle  from  the  stagnation  point.  There  are  statisti¬ 
cal  fluctuations  in  the  DSMC  data  due  to  the  low  vibrational  collision 
probability  at  the  lower  temperatures. 


Vibrational  Temperature  Jump  [K]  Vibrational  Temperature  Jump  [K] 


133 


-  DSMC 

- CFD  (no-slip) 

-  CFD  (Gokcen) 


-  DSMC 

- - CFD  (no-slip) 

-  CFD  (Gokcen) 


(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


-  DSMC 

- CFD  (no-slip) 


-  CFD  (Gokcen) 


-  DSMC 

- - CFD  (no-slip) 

-  CFD  (Gokcen) 


(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  4.61:  Vibrational  temperature  jump  for  a  Mach  25  flow  of  nitrogen  about  a 
cylinder.  $  is  the  angle  from  the  stagnation  point. 
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The  Mach  25,  Kn  =  0.01  and  Kn  =  0.05  cases  also  show  very  reasonable  agreement 
between  the  simulations,  with  the  Lockerby  cases  being  very  close  to  the  DSMC 
solutions.  However,  neither  of  the  CFD  slip  cases  are  able  to  accurately  predict  the 
vibrational  temperature  jump  in  the  Kn  =  0.25  case. 

4.3.7  Computational  Details 

The  computational  details  for  the  nitrogen  simulations  are  given  in  Table  4.8. 

4.4  Summary — Hypersonic  Flow  about  a  Cylinder 

Comparison  of  CFD  and  DSMC  results  for  the  flow  of  argon  about  a  cylinder 
show  that  the  surface  properties  of  pressure,  shear  stress  and  heat  transfer  rates  are 
very  similar  for  the  lower  Knudsen  number  flows  where  the  continuum  hypothesis 
is  valid,  as  expected,  while  the  results  diverge  in  the  higher  Knudsen  number  cases. 
The  surface  pressure  is  least  affected  by  continuum  breakdown,  as  quantified  by  the 
gradient-length  local  Knudsen  number,  among  those  properties  investigated,  and 
seems  to  be  affected  only  by  continuum  breakdown  as  the  shock  and  boundary  layer 
merge  in  the  highest  Knudsen  number  flows.  The  shear  stress  is  most  influenced  by 
nonequilibrium  effects.  The  addition  of  slip  velocity  and  temperature  jump  bound¬ 
ary  conditions  greatly  improve  the  agreement  at  higher  Knudsen  numbers.  Several 
different  types  of  slip  boundary  conditions  are  examined,  and  the  best  agreement 
for  the  surface  properties  is  obtained  when  using  the  generalized  slip  conditions  pro¬ 
posed  by  Gokgen  [29].  With  these  boundary  conditions,  the  differences  in  total  drag 
and  peak  heat  flux  predicted  by  CFD  and  DSMC  increase  from  less  than  1%  at 
Kiioo  =  0.002  to  around  5%  at  Kn^  =  0.25. 

For  the  case  of  a  simple  gas,  the  higher  velocities  associated  with  a  Mach  25 
flow  do  not  seem  to  increase  the  difference  between  the  CFD  and  DSMC  predictions. 
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Table  4.8:  Computational  details  for  a  flow  of  nitrogen  about  a  cylinder.  Total  CPU 
time  is  the  wall  time  multiplied  by  the  number  of  CPUs. 

Mach  10 


Krioo 

DSMC 

Cells 

Particles 

Time  Steps 

CPUs 

Total  CPU  Time* [hours 

0.002 

120,309 

70.9  xlO6 

400,000 

64 

7,296 

0.01 

102,755 

39.4x10s 

300,000 

32 

2,847 

0.05 

35,421 

12.4x10s 

200,000 

16 

460 

0.25 

18,402 

6.2x10s 

200,000 

8 

221 

Klloo 

CFD 

Cells 

Iterations 

CPUs 

Total  CPU  Time*[hours 

0.002 

80,000 

25,000 

8 

300 

0.01 

40,000 

25,000 

8 

150 

0.05 

40,000 

25,000 

8 

150 

0.25 

40,000 

25,000 

8 

150 

Mach  25 


Klloo 

DSMC 

Cells 

Particles 

Time  Steps 

CPUs 

Total  CPU  Time' [hours 

0.002 

102,626 

75.1x10s 

420,000 

64 

8,319 

0.01 

71,997 

37.7x10s 

300,000 

32 

2,734 

0.05 

22,523 

13.1x10s 

200,000 

16 

508 

0.25 

18,395 

6.2x10s 

200,000 

8 

243 

Klloo 

CFD 

Cells 

Iterations 

CPUs 

Total  CPU  Time* [hours 

0.002 

50,000 

45,000 

8 

328 

0.01 

40,000 

25,000 

8 

150 

0.05 

40,000 

25,000 

8 

150 

0.25 

40,000 

25,000 

8 

150 

Approximate 
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Although  the  extent  of  the  region  where  the  continuum  breakdown  parameter  exceeds 
the  critical  value  of  0.05  is  larger  at  the  higher  Mach  number,  the  predicted  surface 
properties  with  the  Gokgen  slip  boundary  conditions  still  remain  well  under  5%  for 
all  but  the  Kn^  =  0.25  case,  where  the  peak  heat  transfer  rates  and  total  drag 
predictions  are  within  6%. 

The  best  agreement  for  the  actual  slip  values  is  obtained  with  the  Lockerby  [49] 
slip  conditions.  While  the  Lockerby  slip  conditions  attempt  to  predict  the  actual 
value  of  the  velocity  slip  and  temperature  jump,  the  Maxwell  slip  conditions  use  a 
higher  velocity  slip  at  the  wall  in  order  to  increase  accuracy  in  flow  properties  further 
from  the  wall.  The  Gokgen  slip  conditions  were  derived  to  accurately  predict  the 
wall  properties  of  shear  stress  and  heat  flux  as  the  Knudsen  number  increases,  rather 
than  the  actual  slip  values. 

Although  there  is  a  significant  amount  of  nonequilibrium  between  the  different 
thermal  modes  (translational,  rotational  and  vibrational)  when  considering  a  flow  of 
nitrogen,  the  trends  are  largely  similar  to  those  noted  when  considering  a  simple  gas 
with  no  internal  degrees  of  freedom.  The  pressure  and  shear  stress  are  least  sensitive 
to  the  nonequilibrium  effects,  while  the  heat  transfer  rate  is  most  sensitive.  Total 
drag  differences  between  CFD  (with  the  best  slip  boundary  conditions)  and  DSMC 
remain  under  3%,  while  peak  heat  flux  differences  are  less  than  8%. 

It  is  also  shown  that  as  the  Knudsen  number  increases,  the  percentage  of  total 
drag  due  to  friction  forces  (versus  pressure)  increases  as  well.  Differences  in  drag 
due  to  skin  friction  also  tend  to  be  larger  than  differences  in  predicted  drag  due 
to  pressure;  thus  the  larger  errors  at  the  higher  Knudsen  numbers  is  due  mostly  to 
errors  in  skin  friction  prediction. 

Differences  in  flow  property  prediction  is  generally  concentrated  in  the  Knudsen 
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layer,  defined  here  as  the  region  of  the  flow  10  mean  free  paths  or  less  from  the  wall 
surface. 


CHAPTER  V 


Hypersonic  Flow  about  a  Wedge 

5.1  Introduction 

The  previous  chapter  considered  the  case  of  a  hypersonic  blunt  body.  This  chapter 
now  considers  the  flow  about  a  wedge  with  a  sharp  leading-edge.  The  cylinder 
simulations  exhibited  a  typical  blunt-body  hypersonic  flow  with  an  unattached  shock, 
followed  by  a  high-temperature,  low  velocity  region  termed  the  shock  layer.  Regions 
of  high  nonequilibrium  were  found  in  the  shock,  the  boundary  layer  and  the  wake.  A 
sharp  leading-edge  body,  on  the  other  hand,  is  characterized  by  an  attached  shock, 
and  supersonic  velocities  throughout  the  flow  (with  the  exception  of  the  boundary 
layer).  Regions  of  high  nonequilibrium  are  expected  near  the  leading-edge  and  in  the 
boundary  layer  and  wake. 

The  wedge  considered  here  has  a  10-deg  half-angle;  the  height  of  the  base  is 
12-inches,  equivalent  to  the  diameter  of  the  cylinder  previously  considered. 

The  inflow  boundary  conditions  are  the  same  as  those  for  the  cylinder;  that  is, 
the  inflow  Mach  numbers  are  10  and  25,  and  the  free  stream  density  of  the  flow  is 
varied  such  that  several  different  regimes  are  considered,  from  the  continuum  through 
the  transitional  to  the  rarefied  regime,  as  shown  in  Table  4.1.  The  global  Knudsen 
number  is  calculated  based  on  the  wedge  base  height,  again  using  the  hard-sphere 
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model  for  the  mean  free-path. 

The  other  relevant  boundary  conditions  are  also  as  shown  in  Table  4.2. 

As  with  the  cylinder,  surface  and  flow  field  properties  for  this  flow  are  presented 
from  two  different  computational  approaches.  First,  CFD  results  are  obtained  using 
the  LeMANS  code.  Previous  simulations  showed  that  the  Gokgen  slip  conditions  [29] 
gives  the  best  agreement  with  DSMC;  therefore,  this  chapter  will  present  only  CFD 
simulations  computed  using  the  Gokgen  boundary  conditions. 

As  with  the  cylinder,  a  mesh-independence  study  is  conducted  in  each  case  to 
determine  the  final  mesh  resolution  used.  For  the  cylinder  it  was  determined  that 
the  wall-normal  spacing  had  the  most  effect  on  mesh-independence  of  the  surface 
properties.  For  the  wedge,  mesh  refined  near  the  leading  edge  in  the  wall-parallel 
direction  also  has  a  large  effect  on  the  surface  properties,  due  to  the  large  gradients 
in  that  area.  For  the  wedge,  in  addition  to  the  wall-normal  mesh  refinement,  the 
number  of  nodes  near  the  leading  edge  is  successively  doubled.  A  mesh-independent 
solution  is  defined  as  one  for  which  the  total  drag  and  peak  heat  transfer  rate  change 
by  1%  or  less  when  using  the  more  refined  mesh. 

Second,  DSMC  results  are  provided  from  the  MONACO  code  for  the  same  flow 
conditions.  Again,  the  mesh  used  for  the  final  solution  is  adapted  such  that  the 
cell  sizes  are  on  the  order  of  one  mean  free  path  or  smaller,  with  the  exceptions  of 
the  Kn  =  0.002  and  Kn  =  0.01  cases  (where  the  subcell  method  is  used  to  select 
particles  for  collisions  to  ensure  physical  accuracy  [4]).  The  flow  near  the  leading- 
edge  exhibits  an  extremely  high  amount  of  nonequilibrium.  Thus,  the  mesh  around 
the  leading-edge  in  each  case  is  adapted  such  that  the  cell  size  is  about  10-40%  of  a 
mean  free  path  in  order  to  sufficiently  resolve  the  flow  details.  Example  meshes  for 
both  DSMC  and  CFD  for  the  different  flow  regimes  are  shown  in  Figure  5.1.  In  the 
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rpsMci 

(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


Figure  5.1:  Example  meshes  for  DSMC  and  CFD  for  the  flow  about  a  wedge.  In  the 
inset,  density  ratio  contours  are  plotted  with  the  mesh,  giving  an  idea  of 
the  mesh  resolution  across  the  shock. 


inset  plots,  density  ratio  contours  are  plotted  with  the  mesh  in  a  detailed  view  of  the 
leading  edge,  giving  an  idea  of  the  mesh  resolution  across  the  shock. 

As  before,  the  DSMC  solutions  are  assumed  to  be  the  more  correct  solutions. 

In  the  results  that  follow,  the  surface  properties  are  presented  in  terms  of  non- 
dimensional  coefficients,  as  defines  in  Eqs.  4.2  -  4.4.  The  surface  properties  in  each 
case  are  plotted  as  a  function  of  the  distance,  S ,  along  the  wedge  surface,  normalized 
by  the  length,  L,  of  the  top  surface.  Thus,  S/L  =  1  is  the  location  of  the  wedge 
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Tabic  5.1:  Total  drag  for  flow  of  argon  about  a  wedge. 


Mach  10 

Mach  25 

Kn^ 

Drag/Length  [N/: 

m]  (%  Difference) 

DSMC 

CFD 

DSMC 

CFD 

0.002 

38.89 

39.42  (1.4%) 

252.3 

257.2  (1.9%) 

0.01 

13.66 

14.22  (4.1%) 

92.74 

99.86  (7.7%) 

0.05 

5.091 

5.639  (10.8%) 

34.39 

40.13  (16.7%) 

0.25 

1.709 

2.051  (20.0%) 

10.35 

13.91  (34.4%) 

shoulder  and  the  beginning  of  the  wake  of  the  flow. 

Along  with  the  surface  properties,  the  maximum  value  of  Kiiqll  at  the  surface 
(based  on  the  CFD  solution)  is  plotted  in  each  case. 

5.2  Argon 

The  flow  of  argon  about  the  wedge  is  considered  first.  Table  5.1  summarizes 
the  total  drag  predicted  by  CFD  and  DSMC.  There  is  reasonable  agrement  at  the 
lowest  Knudsen  number,  with  less  than  2%  difference  for  both  Mach  10  and  Mach 
25.  However,  as  the  flow  becomes  more  rarefied,  the  differences  increase,  with  the 
maximum  differences  of  20%  for  Mach  10  and  over  34%  for  Mach  25  seen  for  the 
Kn  =  0.25  cases. 

Figure  5.2  illustrates  the  percentage  of  total  drag  due  to  pressure  and  friction 
forces,  for  both  DSMC  and  CFD.  It  is  significant  to  note  that  as  the  Knudsen  number 
increases,  the  percentage  of  total  drag  due  to  friction  increases  from  about  50%  at 
Kn  =  0.002  to  over  80%  at  Kn  =  0.25  for  Mach  10.  For  Mach  25,  an  even  larger 
portion  of  the  total  drag  is  due  to  friction  forces — nearly  60%  at  Kn  =  0.002  to  almost 
90%  at  Kn  =  0.25.  This  is  contrasted  with  the  cylinder  in  the  previous  chapter  where 
friction  forces  accounted  for,  at  most,  20%  of  the  total  drag  (at  Kn  =  0.25);  the  vast 
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Figure  5.2:  Percentage  of  total  drag  due  to  pressure  and  friction  for  flow  of  argon 
about  a  wedge.  In  contrast  to  the  cylinder  cases,  here  friction  forces 
account  for  most  of  the  drag. 


majority  of  the  drag  was  due  to  pressure  forces.  It  is  also  significant  to  note,  as 
shown  in  Figure  5.3,  that  the  difference  in  predicted  total  drag  between  CFD  and 
DSMC  is  due  mostly  to  the  differences  predicted  in  the  friction  forces,  as  was  also 
the  case  with  the  cylinder. 

In  all  cases,  CFD  overpredicts  the  total  drag.  Although  these  CFD  simulations 
make  use  of  the  Gokgen  slip  conditions  (which  gave  very  reasonable  results  for  total 
drag  in  the  case  of  the  cylinder),  the  differences  in  total  drag  predictions  here  are 
similar  in  magnitude  to  the  disagreement  produced  when  using  the  no-slip  boundary 
conditions  in  the  case  of  the  cylinder.  This  increased  disagreement  for  the  wedge  is 
simply  due  to  the  fact  that  a  larger  percentage  of  the  total  drag  is  due  to  friction 
forces,  and  the  shear  stress  is  more  sensitive  to  continuum  breakdown  due  to  rarefac¬ 
tion  than  is  pressure;  hence  there  is  more  disagreement  in  the  total  drag  predictions 
for  the  wedge  than  there  was  for  the  cylinder. 

The  peak  heat  transfer  rate  predicted  by  CFD  and  DSMC  is  summarized  in  Table 
5.2,  and  the  differences  are  shown  graphically  in  Figure  5.4.  Here  the  differences 
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(a)  Mach  10  (b)  Mach  25 

Figure  5.3:  Total  drag  difference  from  DSMC  predicted  by  CFD  for  flow  of  argon 
about  a  wedge. 

between  CFD  and  DSMC  for  all  but  the  most  rarefied  cases  are  near  70%.  It  will  be 
shown  below  that  for  the  Kn  =  0.25  cases  the  peak  heating  is  predicted  by  CFD  to 
occur  at  the  shoulder  of  the  wedge  rather  than  at  the  leading  edge  for  these  cases — 
the  differences  in  peak  heating  at  the  leading  edge  is  100%.  In  all  cases,  CFD  also 
underpredicts  the  maximum  heat  flux.  These  differences  are  significantly  larger  than 
were  seen  with  the  cylinder,  and  this  is  the  most  striking  difference  between  surface 
predictions  of  a  blunt  body  compared  with  those  of  a  sharp-leading  edge  body.  It 
will  also  be  shown  below  that  DSMC  predicts  a  very  high  temperature  region  at  the 
leading  edge  in  all  cases;  CFD  is  unable  to  match  these  flow  conditions  due  to  the 
large  effects  of  nonequilibrium  present  near  the  leading  edge. 

5.2.1  Continuum  Breakdown 

As  before,  the  breakdown  parameter  is  calculated  using  both  the  CFD  and  the 
DSMC  solutions  according  to  Equation  3.26.  For  the  case  of  a  wedge  in  a  hypersonic 
flow  of  a  simple  gas,  breakdown  of  the  continuum  hypothesis  is  expected  in  regions  of 


Difference  from  DSMC 


144 


Table  5.2:  Peak  heat  transfer  rate  for  flow  of  argon  about  a  wedge.  The  large  differ¬ 
ences  between  CFD  and  DSMC  are  dne  to  the  failure  of  CFD  to  predict 
the  high  temperatures  at  the  leading  edge. 


Mach  10 

Mach  25 

Klloo 

Peak  Heating  [kW/m2]  (%  Difference) 

DSMC 

CFD 

DSMC 

CFD 

0.002 

239.5 

76.03 

(-68.3%) 

3807. 

1143.  (-70.0%) 

0.01 

47.48 

15.38 

(-67.6%) 

754.2 

231.1  (-69.4%) 

0.05 

9.451 

3.109 

(-67.1%) 

151.0 

46.50  (-69.2%) 

0.25 

1.902 

1.247 

(-100.0%) 

31.76 

26.04  (-100.0%) 

0.002  0.01  0.05  0.25  0.002  0.01  0.05  0.25 


Knudsen  Number 


Knudsen  Number 


(a)  Mach  10 


(b)  Mach  25 


Figure  5.4:  Peak  heat  transfer  rate  difference  from  DSMC  predicted  by  CFD  for  flow 
of  argon  about  a  wedge.  The  large  differences  in  peak  heating  are  due  to 
the  failure  of  CFD  to  predict  the  high  temperatures  at  the  leading  edge. 
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high  gradients  (such  as  the  shock  and  boundary  layer,  and  especially  near  the  leading 
edge)  and  in  regions  of  rarefaction  (such  as  the  wake).  The  amount  of  continuum 
breakdown  is  again  expected  to  increase  as  the  gas  flow  becomes  more  rarefied. 

The  maximum  Kiigll  f°r  each  case  is  plotted  in  Figures  5.5  and  5.6.  The  detail 
of  the  flow  near  the  leading-edge  is  shown  in  the  inset.  In  these  figures,  the  maximum 
gradient  length  local  Knudsen  number  is  computed  from  the  DSMC  (top)  and  CFD 
(bottom)  solutions.  The  light  gray  regions  correspond  to  KnGLL  <  0.05,  dark  gray 
regions  correspond  to  Kiigll  <  1.0,  and  black  regions  correspond  to  Kiigll  >  1.0. 
Notice  that  the  minimum  value  of  Kiigll  f°r  the  black  regions  is  an  order  of  magni¬ 
tude  greater  than  those  shown  for  the  cylinder;  therefore,  the  amount  of  nonequilib¬ 
rium  represented  by  the  black  portions  in  Figures  5.5  and  5.6  is  much  greater  than 
in  those  regions  represented  by  black  in  Figures  4.6  and  4.7. 

For  the  Kn  =  0.002  cases,  there  is  a  significant  amount  of  nonequilibrium  (as 
quantified  by  a  breakdown  parameter  value  exceeding  0.05,  and  represented  by  dark 
gray)  present  in  the  shock,  the  boundary  layer  and  the  wake.  The  size  of  this  region 
of  nonequilibrium  grows  as  the  flow  becomes  rarefied,  until  it  encompasses  nearly  the 
entire  computational  domain  at  Kn  =  0.25.  This  is  very  similar  to  what  was  seen 
in  the  flow  around  the  cylinder.  Note,  however,  the  large  amount  of  nonequilibrium 
near  the  leading-edge,  as  shown  in  the  inset  of  Figures  5.5  and  5.6,  and  near  the 
shoulder  of  the  wedge  as  the  flow  expands  into  the  wake.  As  the  flow  becomes  more 
rarefied,  these  regions  of  severe  nonequilibrium  grow  larger.  It  will  be  seen  below  that 
the  region  of  high  nonequilibrium  near  the  leading-edge  has  a  much  more  significant 
impact  on  the  surface  property  predictions,  even  for  the  cases  that  are  within  the 
continuum  regime,  than  does  the  wake  region. 


146 


rpsMci 

(a)  Kn  =  0.002 


(c)  Kn  =  0.05 


(b)  Kn  =  0.01 


(d)  Kn  =  0.25 


Figure  5.5:  Kiigll  field  for  a  Mach  10  flow  of  argon  about  a  wedge.  The  light  gray 
regions  correspond  to  Kiigll  <  0.05,  dark  gray  regions  correspond  to 
Kiigll  <  TO,  and  black  regions  correspond  to  Kiigll  >  1-0.  Note  that 
the  minimum  value  of  Kiigll  for  the  black  regions  is  an  order  of  magni¬ 
tude  greater  than  those  for  the  cylinder 
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(c)  Kn  =  0.05  (d)  Kn  =  0.25 


Figure  5.6:  Kiigll  field  for  a  Mach  25  flow  of  argon  about  a  wedge.  The  light  gray 
regions  correspond  to  Kiigll  <  0.05,  dark  gray  regions  correspond  to 
Kiigll  <  TO,  and  black  regions  correspond  to  Kiigll  >  1-0.  Note  that 
the  minimum  value  of  Kiigll  for  the  black  regions  is  an  order  of  magni¬ 
tude  greater  than  those  shown  for  the  cylinder. 
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5.2.2  Flow  Field  Properties 

The  density  ratio  fields,  where  the  density  is  normalized  by  the  free  stream  den¬ 
sity,  are  shown  in  Figures  5.7  and  5.8.  The  thin  shock  present  at  Kn  =  0.002  becomes 
more  diffuse  as  the  Knudsen  number  increases  and  there  is  a  low  density  region  in 
the  wake,  causing  breakdown  of  the  continuum  hypothesis.  As  the  Mach  number  in¬ 
creases  from  Mach  10  to  Mach  25,  there  is  little  qualitative  difference  in  the  density 
contours,  although  the  amount  of  compression  through  the  shock  is  higher  at  the 
higher  velocity,  as  shown  by  the  higher  density  ratio  (a  maximum  density  ratio  of 
about  4.5  for  Mach  25  compared  to  about  3.2  for  Mach  10  at  Kn  =  0.002).  DSMC 
also  predicts  a  higher  density  near  the  wedge  surface.  This  is  most  clearly  seen  in 
the  leading-edge  details  and  at  the  higher  Knudsen  numbers. 

The  temperature  fields  predicted  by  both  CFD  and  DSMC  can  be  seen  in  Figures 
5.9  and  5.10.  The  temperature  field  predicted  by  both  methods  seems  to  agree  in  the 
majority  of  the  computational  domain  for  Kn  =  0.002,  Kn  =  0.01  and  even  Kn  =  0.05, 
with  some  differences  in  the  wake.  However,  near  the  leading  edge  (see  insets),  DSMC 
predicts  a  much  higher  temperature  than  CFD  (the  peak  temperatures  are  about 
3,000-3,400  K  for  Mach  10  and  16,000-20,000  K  for  Mach  25).  For  Kn  =  0.002, 
the  DSMC  temperature  at  the  leading  edge  is  about  40%  higher  than  the  CFD 
temperature,  and  the  peak  DSMC  temperature  is  about  15%  higher  than  the  peak 
CFD  temperature.  As  the  Knudsen  number  increases,  the  difference  in  temperature 
at  the  leading  edge  increases  to  nearly  50%  for  Kn  =  0.25.  Although  DSMC  always 
predicts  the  peak  temperature  to  be  at  the  leading  edge,  CFD  predicts  a  higher 
temperature  in  the  wake  for  these  most  rarefied  cases.  It  is  clear  that  CFD  cannot 
accurately  predict  the  temperature  gradients  near  the  leading  edge  even  for  the 
highest  density  cases,  causing  a  large  difference  in  the  predicted  heat  transfer  rates. 
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(c)  Kn  =  0.05  (d)  Kn  =  0.25 


Figure  5.7:  Density  ratio  field  for  a  Mach  10  flow  of  argon  about  a  wedge. 


Mach  25  flow  of  argon  about  a  wedge. 
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(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05  (d)  Kn  =  0.25 


Figure  5.9:  Temperature  field  for  a  Mach  10  flow  of  argon  about  a  wedge.  DSMC 
predicts  a  much  higher  temperature  than  CFD  near  the  leading  edge 
(inset).  Note  that  CFD  predicts  a  higher  temperature  in  the  wake  than 
near  the  leading  edge  for  Kn  =  0.25,  while  DSMC  predicts  a  maximum 
temperature  at  the  leading  edge. 
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(c)  Kn  =  0.05  (d)  Kn  =  0.25 


Figure  5.10:  Temperature  field  for  a  Mach  25  flow  of  argon  about  a  wedge.  DSMC 
predicts  a  much  higher  temperature  than  CFD  near  the  leading  edge 
(inset).  Note  that  CFD  predicts  a  higher  temperature  in  the  wake  than 
near  the  leading  edge  for  Kn  =  0.25,  while  DSMC  predicts  a  maximum 
temperature  at  the  leading  edge. 
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5.2.3  Surface  Properties 

The  surface  property  distributions  (pressure,  shear  stress  and  heat  flux)  for  each 
case  are  now  examined.  The  surface  pressure,  in  the  form  of  a  pressure  coefficient, 
is  shown  in  Figures  5.11  and  5.12.  For  the  blunt-body,  the  surface  pressure  was  the 
least  sensitive  to  nonequilibrium  of  all  the  surface  properties;  for  the  wedge,  however, 
there  are  significant  differences  in  the  pressure  distributions,  even  at  Kn  =  0.002. 
Although  there  are  significant  levels  of  nonequilibrium  in  the  wake  as  well  as  at  the 
leading  edge,  as  shown  by  Figures  5.5  and  5.6,  the  pressure  is  affected  primarily  at 
the  leading  edge.  The  overall  CFD  pressure  distribution  agrees  qualitatively  with 
DSMC  for  all  but  the  Kn  =  0.25  case,  but  the  peak  pressure  at  the  leading  edge  is 
overpredicted  by  CFD.  The  distributions  also  start  to  differ  in  the  wake  for  Kn  =  0.05 
and  Kn  =  0.25,  with  CFD  predicting  a  large  spike  in  the  pressure  at  S/L  =  1  as 
the  flow  begins  to  expand  into  the  wake.  In  all  cases,  CFD  tends  to  overpredict  the 
pressure,  hence  there  is  some  effect  of  the  pressure  on  the  overall  overprediction  of 
total  drag  by  CFD,  although  this  effect  is  not  as  great  as  that  of  shear  stress,  as  was 
shown  in  Figure  5.3. 

The  shear  stress  on  the  wedge  surface  is  seen  in  Figures  5.13  and  5.14.  Here, 
there  is  a  large  spike  in  the  CFD  shear  stress  right  at  the  leading  edge  that  is  not 
shown  in  the  figures.  The  maximum  friction  coefficient  at  that  point  is  annotated  in 
each  case.  This  peak  ranges  from  a  fairly  low  value  of  0.43  for  Mach  25,  Kn  =  0.002 
(which  is  not  much  different  from  the  peak  DSMC  value)  to  a  maximum  of  41.3  for 
Mach  25,  Kn  =  0.25  (compared  to  about  0.36  for  DSMC).  Other  than  this  peak  at 
the  leading  edge,  there  are  fewer  differences  in  the  shear  stress  than  there  were  for 
the  pressure.  However,  the  total  drag  is  affected  most  by  the  friction  forces.  This 
apparent  paradox  is  explained  by  noting  again  the  much  larger  effect  that  shear  stress 
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(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  5.11:  Surface  pressure  coefficient  for  Mach  10  flow  of  argon  about  a  wedge. 

The  maximum  value  of  KncLL  near  the  surface  plotted  on  the  right  axis. 
The  distance  along  the  surface  (including  the  base),  S,  is  normalized  by 
the  top  surface  length,  L. 
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S/L 


(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  5.12:  Surface  pressure  coefficient  for  Mach  25  flow  of  argon  about  a  wedge. 

The  maximum  value  of  KncLL  near  the  surface  plotted  on  the  right  axis. 
The  distance  along  the  surface  (including  the  base),  S,  is  normalized  by 
the  top  surface  length,  L. 
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has  on  the  total  drag  as  seen  in  Figure  5.3. 

The  peak  at  the  leading  edge  does  not  have  any  significant  effect  on  the  total  drag 
due  to  the  small  surface  area  on  which  the  friction  force  acts.  Figures  5.15  and  5.16 
compare  the  accumulated  total  drag  along  the  wedge  surface  due  to  both  friction 
and  pressure  forces.  The  accumulated  total  drag  along  the  surface  is  calculated  by 
summing  the  drag  from  the  leading  edge  up  to  each  point  along  the  surface,  thus 
the  total  drag  is  shown  at  about  S/L  =  1.2.  The  locations  where  there  are  large 
differences  between  the  CFD  and  DSMC  drag  predictions  are  seen  where  the  distance 
between  the  two  lines  on  the  plots  increases.  Thus,  the  differences  in  friction  drag 
for  Kn  =  0.002  and  Kn  =  0.01  occur  along  the  first  20%  of  the  wedge  surface.  For 
Kn  =  0.05,  the  area  where  the  difference  accumulates  most  is  in  the  first  40-50%  of 
the  wedge  surface;  and  for  Kn  =  0.25,  the  differences  accumulate  mostly  between 
20%  and  80%  of  the  wedge  length.  Note  that  the  peak  in  friction  coefficient  at  the 
leading  edge  has  no  significant  impact  on  the  total  drag.  Figures  5.15  and  5.16  also 
demonstrate  that: 

•  The  contribution  of  friction  forces  to  the  total  drag  increases  as  the  density 
decreases. 

•  There  is  no  contribution  to  total  drag  due  to  friction  forces  in  the  wake. 

•  Pressure  forces  on  the  base  of  the  wedge  (in  the  wake)  decrease  the  total  drag. 

•  There  is  little  disagreement  in  the  predictions  of  total  drag  due  to  pressure 
forces  for  most  cases. 

Fairly  significant  disagreement  in  shear  stress  prediction  for  Kn  =  0.05  and 
Kn  =  0.25  is  also  seen  at  S/L  =  1  at  the  sharp  angle  between  the  top  wedge  surface 


and  the  base. 
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(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  5.13:  Surface  friction  coefficient  for  Mach  10  flow  of  argon  about  a  wedge.  The 
maximum  value  of  Kiiqll  near  the  surface  plotted  on  the  right  axis.  The 
distance  along  the  surface  (including  the  base),  S,  is  normalized  by  the 
top  surface  length,  L. 
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(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  5.14:  Surface  friction  coefficient  for  Mach  25  flow  of  argon  about  a  wedge.  The 
maximum  value  of  Kiiqll  near  the  surface  plotted  on  the  right  axis.  The 
distance  along  the  surface  (including  the  base),  S,  is  normalized  by  the 
top  surface  length,  L. 
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(a)  Kn  =  0.002  (b)  Kn  =  0.01 


(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  5.15:  Contributions  of  pressure  and  skin  friction  forces  to  accumulated  total 
drag  for  a  Mach  10  flow  of  argon  about  a  wedge.  The  distance  along  the 
surface  (including  the  base),  S,  is  normalized  by  the  top  surface  length, 
L. 
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(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  5.16:  Contributions  of  pressure  and  skin  friction  forces  to  accumulated  total 
drag  for  a  Mach  25  flow  of  argon  about  a  wedge.  The  distance  along  the 
surface  (including  the  base),  S,  is  normalized  by  the  top  surface  length, 
L. 


161 


The  heat  transfer  rate  distributions  are  shown  in  Figures  5.17  and  5.18.  Previ¬ 
ously,  it  was  shown  that  DSMC  predicts  a  temperature  as  much  as  40%  higher  than 
CFD  at  the  leading  edge.  The  inability  of  CFD  to  model  the  flow  at  the  leading  edge 
has  a  large  effect  on  the  heat  flux.  Although  the  qualitative  agreement  is  good  at 
Kn  =  0.002,  DSMC  predicts  a  heating  coefficient  of  nearly  0.19,  while  CFD  predicts 
a  much  lower  heating  coefficient  of  about  0.07.  As  the  Knudsen  number  increases, 
the  DSMC  peak  heating  coefficient  remains  at  about  0.19  -  0.20,  while  the  CFD  peak 
heating  coefficient  remains  near  0.05  -  0.06.  Thus,  the  difference  in  heating  rate  is 
around  70%  for  most  of  the  cases.  For  Kn  =  0.25,  CFD  predicts  a  peak  heating  rate 
at  the  shoulder,  with  no  heating  predicted  at  the  leading  edge.  The  actual  error  in 
heating  rate  prediction,  then,  is  much  greater  than  the  35%  and  18%  cited  in  Table 
5.2  and  shown  in  Figure  5.4. 

5.2.4  Slip  Quantities 

The  velocity  slip  along  the  wedge  surface  is  shown  in  Figures  5.19  and  5.20.  In  all 
cases,  CFD  predicts  a  peak  velocity  slip  of  about  2000  m/s  (for  Mach  10)  and  about 
5000  m/s  (for  Mach  25)  near  the  leading  edge,  while  DSMC  predicts  a  maximum 
velocity  slip  of  about  700  m/s  (for  Mach  10)  and  about  1100  m/s  (for  Mach  25). 
Past  the  leading  edge,  the  velocity  slip  is  very  quickly  reduced  to  a  nearly  constant 
finite  value  until  the  wedge  shoulder,  where  a  sharp  increase  is  seen.  For  Kn  =  0.002 
Kn  =  0.01  and  even  Kn  =  0.05  to  some  extent,  CFD  qualitatively  agrees  fairly  well 
with  DSMC.  However,  for  Kn  =  0.25  this  agreement  worsens  considerably.  The 
locations  where  the  disagreement  is  most  apparent,  at  and  near  the  leading  edge,  are 
also  the  locations  where  the  shear  stress  distributions  differ  the  most.  Thus,  there  is 
a  good  correlation  between  velocity  slip  disagreement  and  shear  stress  disagreement. 
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S/L 


(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


S/L 


S/L 


(c)  Kn  =  0.05  (d)  Kn  =  0.25 

Figure  5.17:  Surface  heating  coefficient  for  Mach  10  flow  of  argon  about  a  wedge.  The 
maximum  value  of  Kiiqll  near  the  surface  plotted  on  the  right  axis.  The 
distance  along  the  surface  (including  the  base),  S,  is  normalized  by  the 
top  surface  length,  L. 
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(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


S/L 


S/L 


(c)  Kn  =  0.05  (d)  Kn  =  0.25 

Figure  5.18:  Surface  heating  coefficient  for  Mach  25  flow  of  argon  about  a  wedge.  The 
maximum  value  of  Kiiqll  near  the  surface  plotted  on  the  right  axis.  The 
distance  along  the  surface  (including  the  base),  S,  is  normalized  by  the 
top  surface  length,  L. 
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DSMC 


(a)  Kn  =  0.002 


(c)  Kn  =  0.05 


DSMC 


(b)  Kn  =  0.01 


(d)  Kn  =  0.25 


Figure  5.19:  Velocity  slip  for  a  Mach  10  flow  of  argon  about  a  wedge.  The  distance 
along  the  surface  (including  the  base),  S,  is  normalized  by  the  top 
surface  length,  L. 
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DSMC 


(a)  Kn  =  0.002 


DSMC 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05  (d)  Kn  =  0.25 


Figure  5.20:  Velocity  slip  for  a  Mach  25  flow  of  argon  about  a  wedge.  The  distance 
along  the  surface  (including  the  base),  S,  is  normalized  by  the  top 
surface  length,  L. 
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The  temperature  jump  is  shown  in  Figures  5.21  and  5.22.  The  DSMC  profile  of 
temperature  jump  is  very  similar  to  that  for  velocity  slip;  the  temperature  jump  is 
highest  at  the  leading  edge,  with  a  gradual  decrease  to  an  almost  constant  value  and 
a  slight  bump  around  the  wedge  shoulder.  The  peak  temperature  jump  predicted 
by  DSMC  is  between  about  2000-2200  K  (for  Mach  10)  and  10000-11000  K  (for 
Mach  25).  CFD  agrees  fairly  well  for  Kn  =  0.002,  but  shows  large  disagreements  for 
all  other  cases,  particularly  in  the  wake.  The  CFD  peak  temperature  jump  is  not 
much  different  from  that  of  DSMC,  especially  for  the  Mach  25  cases.  There  does  not 
seem  to  be  any  strong  correlation  between  the  temperature  jump  agreement  and  the 
disagreement  between  heat  transfer  rates. 

5.2.5  Computational  Details 

The  computational  details  for  the  simulations  of  a  hypersonic  flow  of  argon  about 
a  wedge  discussed  in  this  section  are  shown  in  Table  5.3. 
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DSMC 


S/L 


(a)  Kn  =  0.002 


(c)  Kn  =  0.05 


DSMC 


(b)  Kn  =  0.01 


(d)  Kn  =  0.25 


Figure  5.21:  Temperature  jump  for  a  Mach  10  flow  of  argon  about  a  wedge.  The 
distance  along  the  surface  (including  the  base),  S,  is  normalized  by  the 
top  surface  length,  L. 
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DSMC 


S/L 


(a)  Kn  =  0.002 


DSMC 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05  (d)  Kn  =  0.25 


Figure  5.22:  Temperature  jump  for  a  Mach  25  flow  of  argon  about  a  wedge.  The 
distance  along  the  surface  (including  the  base),  S,  is  normalized  by  the 
top  surface  length,  L. 


169 


Table  5.3:  Computational  details  for  a  flow  of  argon  about  a  wedge.  Total  CPU  time 
is  the  wall  time  multiplied  by  the  number  of  CPUs. 

Mach  10 


Klloo 

DSMC 

Cells 

Particles 

Time  Steps 

CPUs 

Total  CPU  Time* [hours 

0.002 

187,036 

50.7x10s 

260,000 

32 

3,072 

0.01 

155,126 

47.6  x10s 

200,000 

32 

2,220 

0.05 

26,329 

9.6x10s 

200,000 

8 

636 

0.25 

8,083 

5.3x10s 

200,000 

4 

159 

Klloo 

CFD 

Cells 

Iterations 

CPUs 

Total  CPU  Time*[hours 

0.002 

25,200 

30,000 

8 

56 

0.01 

25,200 

20,000 

8 

38 

0.05 

28,400 

30,000 

8 

64 

0.25 

30,175 

30,000 

8 

66 

Mach  25 


Klloo 

DSMC 

Cells 

Particles 

Time  Steps 

CPUs 

Total  CPU  Time' [hours 

0.002 

107,421 

32.4x10s 

360,000 

32 

3,494 

0.01 

94,596 

24.2  x10s 

200,000 

16 

1,275 

0.05 

19,398 

4.2x10s 

200,000 

4 

125 

0.25 

6,110 

5.3x10s 

200,000 

4 
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Klloo 

CFD 

Cells 

Iterations 

CPUs 

Total  CPU  Time*[hours 

0.002 

29,050 

30,000 

8 

60 

0.01 

25,200 

30,000 

8 

66 

0.05 

28,400 

30,000 

8 

60 

0.25 

30,175 

30,000 

8 

80 

Approximate 
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Tabic  5.4:  Total  drag  for  flow  of  nitrogen  about  a  wedge. 


Mach  10 

Mach  25 

Klloo 

Drag/Length  [N/: 

m]  (%  Difference) 

DSMC 

CFD 

DSMC 

CFD 

0.002 

29.92 

30.34 

(1.4%) 

183.8 

189.9  (3.3%) 

0.01 

10.48 

10.86 

(3.6%) 

67.82 

72.34  (6.7%) 

0.05 

3.946 

4.255 

(7.8%) 

26.10 

29.03  (11.2%) 

0.25 

1.364 

1.547 

(13.4%) 

8.488 

10.05  (18.4%) 

5.3  Nitrogen 

The  flow  of  nitrogen  about  the  wedge  is  now  considered.  Table  5.4  summarizes 
the  total  drag  predicted  by  CFD  and  DSMC.  Again,  there  is  reasonable  agrement  at 
the  lowest  Knudsen  number,  with  less  than  2%  difference  for  Mach  10,  but  with  a 
little  more  than  3%  for  Mach  25.  As  the  flow  becomes  more  rarefied,  the  differences 
increase,  although  the  maximum  differences  of  about  13%  and  18%  at  Kn  =  0.25  are 
lower  than  that  seen  with  argon. 

Figure  5.23  again  illustrates  the  percentage  of  total  drag  due  to  pressure  and 
friction  forces,  for  both  DSMC  and  CFD.  The  amount  of  drag  due  to  friction  again 
ranges  from  around  50%  for  Mach  10,  Kn  =  0.002,  to  about  90%  for  Mach  25, 
Kn  =  0.25.  A  slightly  larger  percentage  of  the  drag  is  due  to  friction  at  the  higher 
velocity,  as  was  seen  with  argon,  and  the  contribution  of  friction  forces  to  the  total 
drag  also  increases  with  increasing  Knudsen  number.  It  is  also  shown  in  Figure 
5.24,  that  the  difference  in  predicted  total  drag  between  CFD  and  DSMC  is  due 
mostly  to  the  differences  predicted  in  the  friction  forces.  Note,  however,  that  the 
differences  due  to  pressure  drag  are  much  lower  than  was  the  case  with  argon;  for 
instance,  the  argon  case  for  Mach  25,  Kn  =  0.25  case  had  a  difference  of  over  5% 
in  the  pressure  drag  prediction  (see  Figure  5.3),  while  the  currently  considered  cases 
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Figure  5.23:  Percentage  of  total  drag  due  to  pressure  and  friction  for  flow  of  nitrogen 
about  a  wedge.  In  contrast  to  the  cylinder  cases,  here  friction  forces 
account  for  most  of  the  drag. 


show  a  maximum  of  only  about  2%  for  the  same  Mach  and  Knudsen  numbers.  The 
pressure  drag  differences  for  the  Mach  10  nitrogen  cases  are  insignificant. 

Peak  heating  rates  for  the  nitrogen  cases  are  summarized  in  Table  5.5,  and  il¬ 
lustrated  graphically  in  Figure  5.25.  Similar  to  the  argon  cases,  the  error  in  peak 
heating  is  about  70%  for  most  cases,  and  decreases  to  about  40%  and  25%  for  Mach 
10  and  Mach  25,  respectively.  Again,  though,  the  decrease  in  peak  heating  differ¬ 
ences  is  due  to  CFD’s  prediction  of  peak  heating  occurring  at  the  wedge  shoulder 
rather  than  at  the  leading  edge — the  actual  error  at  the  leading  edge  is  nearly  100%. 

5.3.1  Continuum  Breakdown 

Again,  the  breakdown  parameter  is  calculated  using  both  the  CFD  and  the  DSMC 
solutions  according  to  Equation  3.26.  Continuum  breakdown  is  expected  near  the 
leading  edge  and  in  the  wake,  as  was  the  case  for  the  flow  of  argon,  with  the  degree 
of  nonequilibrium  increasing  with  increasing  Knudsen  number. 

The  maximum  Kiigll  for  each  case  is  plotted  in  Figures  5.26  and  5.27,  with 


Difference  from  DSMC 
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Knudsen  Number  Knudsen  Number 

(a)  Mach  10  (b)  Mach  25 

Figure  5.24:  Total  drag  difference  from  DSMC  predicted  by  CFD  for  flow  of  nitrogen 
about  a  wedge. 


Table  5.5:  Peak  heat  transfer  rate  for  flow  of  nitrogen  about  a  wedge.  The  large 
differences  between  CFD  and  DSMC  are  due  to  the  failure  of  CFD  to 
predict  the  high  temperatures  at  the  leading  edge. 


Mach  10 

Mach  25 

Klloo 

Peak  Heating  [kW / nr2]  (%  Difference) 

DSMC 

CFD 

DSMC 

CFD 

0.002 

209.7 

56.76  (-72.9%) 

3333. 

869.5  (-73.9%) 

0.01 

41.84 

11.46  (-72.6%) 

663.1 

172.9  (-73.9%) 

0.05 

8.329 

2.308  (-72.3%) 

132.8 

36.62  (-72.4%) 

0.25 

1.694 

1.032  (-39.1%) 

27.31 

20.89  (-23.5%) 
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Knudsen  Number 


Knudsen  Number 


(a)  Mach  10 


(b)  Mach  25 


Figure  5.25:  Peak  heat  transfer  rate  difference  from  DSMC  predicted  by  CFD  for 
flow  of  nitrogen  about  a  wedge.  The  large  differences  in  peak  heating 
are  due  to  the  failure  of  CFD  to  predict  the  high  temperatures  at  the 
leading  edge. 


the  detail  of  the  leading  edge  in  the  inset  of  each  frame.  The  maximum  gradient 
length  local  Knudsen  number  is  computed  from  the  DSMC  (top)  and  CFD  (bottom) 
solutions.  As  before,  the  light  gray  regions  correspond  to  Kiigll  <  0.05,  dark  gray 
regions  correspond  to  Kiiqll  <  1-0  and  black  regions  correspond  to  KiiGll  >  1-0. 

The  locations  and  values  of  the  breakdown  parameter  are  very  similar  to  what  was 
seen  for  argon.  A  large  degree  of  nonequilibrium  is  predicted  near  the  leading  edge 
and  the  wedge  shoulder,  and  as  the  flow  becomes  more  rarefied,  the  regions  where 
the  continuum  breakdown  parameter  exceeds  the  critical  value  of  0.05  grow  until 
nearly  all  of  the  computational  domain  is  expected  to  require  DSMC  for  accurate 
modeling. 

5.3.2  Flow  Field  Properties 

The  density  ratio  fields,  where  the  density  is  normalized  by  the  free  stream  den¬ 
sity,  are  shown  in  Figures  5.28  and  5.29.  The  density  ratio  results  are  very  similar 
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(c)  Kn  =  0.05  (d)  Kn  =  0.25 


Figure  5.26:  Kiiqll  field  for  a  Mach  10  flow  of  nitrogen  about  a  wedge.  The  light 
gray  regions  correspond  to  Kiigll  <  0.05,  dark  gray  regions  correspond 
to  Kiigll  <  TO,  and  black  regions  correspond  to  Kiigll  >  TO.  Note 
that  the  minimum  value  of  Kiigll  for  the  black  regions  is  an  order  of 
magnitude  greater  than  those  for  the  cylinder. 
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rpsMci 

(a)  Kn  =  0.002 


(c)  Kn  =  0.05 


fpslvicl 

(b)  Kn  =  0.01 


(d)  Kn  =  0.25 


Figure  5.27:  Kiiqll  field  for  a  Mach  25  flow  of  nitrogen  about  a  wedge.  The  light 
gray  regions  correspond  to  Kiigll  <  0.05,  dark  gray  regions  correspond 
to  Kiigll  <  TO,  and  black  regions  correspond  to  Kiigll  >  TO.  Note 
that  the  minimum  value  of  Kiigll  for  the  black  regions  is  an  order  of 
magnitude  greater  than  those  for  the  cylinder. 
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(a)  Kn  =  0.002 


(c)  Kn  =  0.05 


(b)  Kn  =  0.01 


(d)  Kn  =  0.25 


Figure  5.28:  Density  ratio  field  for  a  Mach  10  flow  of  nitrogen  about  a  wedge. 


to  those  for  argon,  with  the  exception  that  nitrogen  achieves  a  higher  compression 
ratio  through  the  shock;  the  maximum  density  ratio  is  about  3.8  for  Mach  10  and 
about  5  for  Mach  25.  The  shock  here  is  not  normal  (as  with  the  cylinder),  and  the 
compression  ratio  is  not  as  high  as  theoretically  possible;  nevertheless,  a  diatomic 
gas  can  be  compressed  more  than  a  monatomic  gas.  DSMC  again  predicts  a  higher 
density  near  the  wedge  surface,  as  seen  in  the  leading-edge  details  and  at  the  higher 
Knudsen  numbers. 

The  translational/rotational  temperature  fields  predicted  by  both  CFD  and  DSMC 


(c)  Kn  =  0.05 


Figure  5.29:  Density  ratio  field 
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can  be  seen  in  Figures  5.30  and  5.31.  Again,  the  area  of  most  concern  is  near  the 
leading  edge  where  DSMC  predicts  a  higher  temperature  than  does  CFD.  Although 
the  peak  temperatures  in  the  nitrogen  flows  are  lower  than  those  in  the  argon  flows 
(around  1,700-1,900  K  for  Mach  10  and  9,000-10,000  K  for  Mach  25),  CFD  again 
underpredicts  the  temperature  by  about  40-50%.  Similarly  to  argon,  DSMC  always 
predicts  the  peak  temperature  to  be  at  the  leading  edge,  while  CFD  predicts  a  higher 
temperature  in  the  wake  for  these  most  rarefied  cases. 

Due  to  the  additional  internal  energy  modes  present  in  a  diatomic  gas  such  as 
nitrogen,  a  vibrational  temperature  is  also  modeled  for  the  wedge.  The  vibrational 
temperature  fields  are  shown  in  Figures  5.32  and  5.33.  For  Mach  10,  the  maximum 
vibrational  temperature  is  near  the  wall  temperature  of  500  K,  again  indicating 
(along  with  the  relatively  low  translational/rotational  temperatures  of  approximately 
2,000  K)  that  the  vibrational  modes  are  only  activated  due  to  the  wall  boundary 
conditions.  The  effects  of  the  CFD  temperature  jump  conditions  are  seen  in  Figure 
5.32  as  the  vibrational  temperature  does  not  attain  the  wall  temperature  of  500  K 
right  at  the  leading  edge. 

The  Mach  25  cases  do  exhibit  some  vibrational  excitation  due  to  elevated  tem¬ 
peratures;  the  peak  vibrational  temperature  is  about  1,800  K  for  the  Mach  25, 
Kn  =  0.002  case.  Here  it  is  also  seen  that  DSMC  predicts  a  higher  level  of  vibrational 
activation  than  does  CFD.  This  is  most  likely  due  to  the  method  used  to  adjust  the 
DSMC  vibrational  collision  probability  to  match  the  theoretical  vibrational  collision 
probability,  as  explained  in  Chapter  Ill.  Recall  that  the  DSMC  vibrational  collision 
probability  is  multiplied  by  a  correction  factor  based  on  the  maximum  translational 
temperature  seen  in  the  flow  field.  Here,  the  maximum  translational  temperatures 
are  approximately  10,000  K,  and  a  correction  factor  of  1.43  is  used  (see  Table  3.2). 
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(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  5.30:  Translational/rotational  temperature  field  for  a  Mach  10  flow  of  nitro¬ 
gen  about  a  wedge.  DSMC  predicts  a  much  higher  temperature  than 
CFD  near  the  leading  edge  (inset).  Note  that  CFD  predicts  a  higher 
temperature  in  the  wake  than  near  the  leading  edge  for  Kn  =  0.25,  while 
DSMC  predicts  a  maximum  temperature  at  the  leading  edge. 
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rpsMci 

(a)  Kn  =  0.002 


rpslvicl 

(b)  Kn  =  0.01 


(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  5.31:  Translational/rotational  temperature  field  for  a  Mach  25  flow  of  nitro¬ 
gen  about  a  wedge.  DSMC  predicts  a  much  higher  temperature  than 
CFD  near  the  leading  edge  (inset).  Note  that  CFD  predicts  a  higher 
temperature  in  the  wake  than  near  the  leading  edge  for  Kn  =  0.25,  while 
DSMC  predicts  a  maximum  temperature  at  the  leading  edge. 
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(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  5.32:  Vibrational  temperature  field  for  a  Mach  10  flow  of  nitrogen  about  a 
wedge.  Note  that  the  vibrational  activation  is  due  to  the  wall  vibrational 
temperature  of  500  K. 


However,  the  maximum  vibrational  temperatures  are  not  near  the  leading  edge  of 
the  wedge,  but  further  back  along  the  surface  where  the  translational/rotational 
temperatures  are  much  lower,  around  4,000-5,000  K.  The  DSMC  vibrational  colli¬ 
sion  probability  for  temperatures  in  that  range  are  higher  than  the  theoretical  model 
(again  see  Table  3.2  and  Figure  3.2);  thus  the  correction  factor  should  be  smaller 
than  unity  (around  0.79),  rather  than  larger.  Nevertheless,  vibrational  temperature 
differences  do  not  seem  to  affect  the  surface  properties  significantly. 
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(b)  Kn  =  0.01 


(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  5.33:  Vibrational  temperature  field  for  a  Mach  25  flow  of  nitrogen  about  a 
wedge.  Note  that  the  wall  vibrational  temperature  is  1500  K. 
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5.3.3  Surface  Properties 

The  surface  property  distributions  (pressure,  shear  stress  and  heat  flux)  for  each 
of  the  cases  are  now  examined.  The  surface  pressure,  in  the  form  of  a  pressure 
coefficient,  is  shown  in  Figures  5.34  and  5.35.  The  overall  CFD  pressure  distribution 
agrees  qualitatively  with  DSMC  for  all  but  the  Kn  =  0.25  case,  but  the  peak  pressure 
at  the  leading  edge  is  overpredicted  by  CFD.  The  distributions  also  start  to  differ 
in  the  wake  for  Kn  =  0.05  and  Kn  =  0.25,  with  CFD  predicting  a  large  spike  in 
the  pressure  at  S/L  =  1  as  the  flow  begins  to  expand  into  the  wake.  In  all  cases, 
CFD  tends  to  overpredict  the  pressure  and  there  is  some  effect  of  the  pressure  on 
the  overall  overprediction  of  total  drag  by  CFD,  although  this  effect  is  not  as  great 
as  that  of  shear  stress,  as  was  shown  in  Figure  5.24. 

The  shear  stress  on  the  wedge  surface  is  seen  in  Figures  5.36  and  5.37.  The  large 
spike  in  the  CFD  shear  stress  at  the  leading  edge  does  not  significantly  affect  the 
total  drag,  and  the  value  of  this  peak  is  annotated  on  the  plots.  Other  than  this 
peak  at  the  leading  edge,  there  are  fewer  differences  in  the  shear  stress  than  there 
were  for  the  pressure.  However,  the  total  drag  is  affected  most  by  the  friction  forces. 
Again,  this  is  explained  by  noting  the  much  larger  effect  that  shear  stress  has  on  the 
total  drag  as  seen  in  Figure  5.3. 

Figures  5.38  and  5.39  compare  the  accumulated  total  drag  along  the  wedge  surface 
due  to  both  friction  and  pressure  forces  for  the  nitrogen  cases.  Thus,  the  differences 
in  friction  drag  for  Kn  =  0.002  and  Kn  =  0.01  occur  along  the  first  20%  of  the  wedge 
surface.  For  Kn  =  0.05,  the  area  where  the  difference  accumulates  most  is  in  the  first 
40-50%  of  the  wedge  surface;  and  for  Kn  =  0.25,  the  differences  accumulate  mostly 
between  20%  and  80%  of  the  wedge  length.  The  same  trends  are  noted  here  with 
nitrogen  as  they  were  with  argon: 


Pressure  Coefficient,  Cp  Pressure  Coefficient,  C 
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(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  5.34:  Surface  pressure  coefficient  for  a  Mach  10  flow  of  nitrogen  about  a 
wedge.  The  maximum  value  of  KnQLL  near  the  surface  plotted  on  the 
right  axis.  The  distance  along  the  surface  (including  the  base),  S,  is 
normalized  by  the  top  surface  length,  L. 
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S/L 


(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  5.35:  Surface  pressure  coefficient  for  a  Mach  25  flow  of  nitrogen  about  a 
wedge.  The  maximum  value  of  KnQLL  near  the  surface  plotted  on  the 
right  axis.  The  distance  along  the  surface  (including  the  base),  S,  is 
normalized  by  the  top  surface  length,  L. 
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(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  5.36:  Surface  friction  coefficient  for  a  Mach  10  flow  of  nitrogen  about  a  wedge. 

The  maximum  value  of  Kiigll  near  the  surface  plotted  on  the  right  axis. 
The  distance  along  the  surface  (including  the  base),  S,  is  normalized  by 
the  top  surface  length,  L. 
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(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  5.37:  Surface  friction  coefficient  for  a  Mach  25  flow  of  nitrogen  about  a  wedge. 

The  maximum  value  of  Kiigll  near  the  surface  plotted  on  the  right  axis. 
The  distance  along  the  surface  (including  the  base),  S,  is  normalized  by 
the  top  surface  length,  L. 
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•  The  large  value  of  skin  friction  coefficient  at  the  leading  edge  does  not  signifi¬ 
cantly  affect  the  total  drag. 

•  The  contribution  of  friction  forces  to  the  total  drag  increases  as  the  density 
decreases. 

•  There  is  no  contribution  to  total  drag  due  to  friction  forces  in  the  wake. 

•  Pressure  forces  on  the  base  of  the  wedge  (in  the  wake)  decrease  the  total  drag. 

•  There  is  little  disagreement  in  the  predictions  of  total  drag  due  to  pressure 
forces  for  most  cases. 

Additionally,  the  good  agreement  between  drag  due  to  pressure,  as  Figure  5.24  shows, 
is  not  due  to  perfect  agreement  between  the  pressure  distributions  along  the  entire 
surface  (as  discussed  above).  Instead,  the  accumulated  drag  clue  to  pressure  does 
show  some  differences  along  the  surface,  but  the  final  total  values  happen  to  agree. 

The  heat  transfer  rate  distributions  are  shown  in  Figures  5.40  and  5.41.  Previ¬ 
ously,  it  was  shown  that  DSMC  predicts  a  much  higher  temperature  at  the  leading 
edge.  The  inability  of  CFD  to  model  the  flow  at  the  leading  edge  again  has  a  large 
effect  on  the  heat  flux.  The  same  general  trends  with  argon  are  seen  here,  although 
a  larger  heat  flux  coefficient  is  predicted  for  nitrogen.  The  peak  DSMC  heating  co¬ 
efficients  remain  near  0.22,  and  the  peak  CFD  heating  coefficients  remain  well  under 
0.1.  Thus,  the  difference  in  heating  rate  is  also  around  70%  for  most  of  the  cases. 
For  Kn  =  0.25,  CFD  predicts  a  peak  heating  rate  at  the  shoulder,  with  no  heating 
predicted  at  the  leading  edge.  Thus,  the  actual  error  in  heating  rate  prediction  is 
much  greater  than  the  40%  and  25%  cited  in  Table  5.5  and  shown  in  Figure  5.25. 
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(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  5.38:  Contributions  of  pressure  and  skin  friction  forces  to  accumulated  total 
drag  for  a  Mach  10  flow  of  nitrogen  about  a  wedge.  The  distance  along 
the  surface  (including  the  base),  S,  is  normalized  by  the  top  surface 
length,  L. 
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(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  5.39:  Contributions  of  pressure  and  skin  friction  forces  to  accumulated  total 
drag  for  a  Mach  25  flow  of  nitrogen  about  a  wedge.  The  distance  along 
the  surface  (including  the  base),  S,  is  normalized  by  the  top  surface 
length,  L. 
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(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05 


(d)  Kn  =  0.25 


Figure  5.40:  Surface  heating  coefficient  for  a  Mach  10  flow  of  nitrogen  about  a  wedge. 

The  maximum  value  of  KncLL  near  the  surface  plotted  on  the  right  axis. 
The  distance  along  the  surface  (including  the  base),  S,  is  normalized  by 
the  top  surface  length,  L. 
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S/L 


(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05  (d)  Kn  =  0.25 

Figure  5.41:  Surface  heating  coefficient  for  a  Mach  25  flow  of  nitrogen  about  a  wedge. 

The  maximum  value  of  KncLL  near  the  surface  plotted  on  the  right  axis. 
The  distance  along  the  surface  (including  the  base),  S,  is  normalized  by 
the  top  surface  length,  L. 
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5.3.4  Slip  Quantities 

The  velocity  slip  along  the  wedge  surface  is  shown  in  Figures  5.42  and  5.43. 
Overall,  the  velocity  slip  profiles  are  very  similar  to  those  obtained  with  argon.  Here, 
CFD  predicts  a  peak  velocity  slip  between  2,000-2,500  m/s  for  Mach  10  and  between 
5000-6,000  m/s  for  Mach  25  near  the  leading  edge,  while  DSMC  predicts  a  maximum 
velocity  slip  of  about  750  m/s  for  Mach  10  and  about  1,500  m/s  for  Mach  25.  For 
Kn  =  0.002,  Kn  =  0.01  and  even  Kn  =  0.05  to  some  extent,  CFD  qualitatively  agrees 
fairly  well  with  DSMC.  However,  for  Kn  =  0.25  this  agreement  worsens  considerably. 
Once  again,  the  locations  where  the  disagreement  is  most  apparent,  at  and  near  the 
leading  edge,  are  also  the  locations  where  the  shear  stress  distributions  differ  the 
most.  Thus,  there  is  a  good  correlation  between  velocity  slip  disagreement  and  shear 
stress  disagreement. 

The  translational/rotational  temperature  jump  profiles  for  the  nitrogen  flows, 
shown  in  Figures  5.44  and  5.45,  are  qualitatively  similar  to  the  temperature  jump 
profiles  obtained  for  argon.  Also  included  here  are  the  translational  and  rotational 
temperature  jump  values  from  the  DSMC  simulations.  Although  the  peak  transla¬ 
tional  temperature  jump  values  predicted  by  DSMC  are  nearly  2,000  K  for  Mach  10 
and  nearly  10,000  K  for  Mach  25,  which  are  near  those  predicted  for  argon,  there 
is  significant  thermal  nonequilibrium  and  the  rotational  temperature  jump  is  much 
lower  than  the  translational  temperature  jump;  at  the  leading  edge,  the  rotational 
temperature  jump  is  negative.  For  Kn  =  0.05  and  Kn  =  0.25,  the  translational  and 
rotational  temperatures  along  the  surface  never  do  fully  equilibrate,  although  they 
are  closer  to  equilibrium  at  the  base  of  the  wedge.  Nevertheless,  CFD  agrees  moder¬ 
ately  well  with  the  averaged  translational/rotational  temperature  jump  predicted  by 
DSMC  for  the  lower  Knudsen  number  flows  for  Mach  10;  there  is  more  disagreement 
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DSMC 


(a)  Kn  =  0.002 


DSMC 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05  (d)  Kn  =  0.25 


Figure  5.42:  Velocity  slip  for  a  Mach  10  flow  of  nitrogen  about  a  wedge.  The  distance 
along  the  surface  (including  the  base),  S,  is  normalized  by  the  top 
surface  length,  L. 
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(c)  Kn  =  0.05 
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Figure  5.43:  Velocity  slip  for  a  Mach  25  flow  of  nitrogen  about  a  wedge.  The  distance 
along  the  surface  (including  the  base),  S,  is  normalized  by  the  top 
surface  length,  L. 
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for  Mach  25. 

The  vibrational  temperature  jump  is  also  seen  in  Figures  5.46  and  5.47.  For 
Mach  10,  there  is  very  little  vibrational  excitation  and  the  vibrational  temperature 
jump  predicted  by  both  methods  is  near  zero,  except  at  the  leading  edge  where  CFD 
overpredicts  the  amount  of  vibrational  temperature  slip.  The  Kn  =  0.25  case  is 
an  exception,  where  DSMC  predicts  a  small,  negative,  temperature  jump  value.  For 
Mach  25,  the  situation  is  very  similar  to  that  for  Mach  10,  except  that  the  differences 
between  CFD  and  DSMC  are  more  pronounced  for  Kn  =  0.01  and  Kn  =  0.05. 

5.3.5  Computational  Details 

The  computational  details  for  the  simulations  of  a  hypersonic  nitrogen  flow  about 
a  wedge  discussed  in  this  section  are  shown  in  Table  5.6. 

5.4  Summary — Hypersonic  Flow  about  a  Wedge 

The  sharp-leading  edge  geometry  of  the  wedge  leads  to  additional  flow  phenomena 
not  seen  with  the  cylinder,  which  in  turn  affect  the  surface  property  predictions. 

The  differences  in  total  drag  predicted  by  CFD  and  DSMC  are  greater  than  they 
were  with  the  cylinder  geometry,  due  to  the  larger  effect  of  the  friction  forces  on  the 
drag.  Additionally,  there  are  greater  differences  in  the  pressure  profiles  along  the 
wedge  surface,  but  the  effect  on  total  drag  is  relatively  small  (due  to  the  small  angle 
of  the  surface  with  the  flow). 

The  amount  of  nonequilibrium  near  the  leading-edge  significantly  affects  the  pre¬ 
diction  of  temperature  gradients,  and  thus  there  are  significant  differences  in  the 
heat  transfer  rate  predictions-CFD  fails  to  adequately  predict  the  large  heat  fluxes 
near  the  leading  edge  of  an  inhnitely-sharp  wedge. 
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DSMC 


DSMC 


(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05  (d)  Kn  =  0.25 

Figure  5.44:  Translational/rotational  temperature  jump  for  a  Mach  10  flow  of  nitro¬ 
gen  about  a  wedge.  The  separate  translational  and  rotational  tempera¬ 
ture  jumps  from  DSMC  are  plotted  along  with  the  average  temperature 
jump  from  CFD  and  DSMC.  The  distance  along  the  surface  (including 
the  base),  S,  is  normalized  by  the  top  surface  length,  L. 
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DSMC 


(a)  Kn  =  0.002 


(b)  Kn  =  0.01 


(c)  Kn  =  0.05  (d)  Kn  =  0.25 

Figure  5.45:  Translational/rotational  temperature  jump  for  a  Mach  25  flow  of  nitro¬ 
gen  about  a  wedge.  The  separate  translational  and  rotational  tempera¬ 
ture  jumps  from  DSMC  are  plotted  along  with  the  average  temperature 
jump  from  CFD  and  DSMC.  The  distance  along  the  surface  (including 
the  base),  S,  is  normalized  by  the  top  surface  length,  L. 
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Figure  5.46:  Vibrational  temperature  jump  for  a  Mach  10  flow  of  nitrogen  about 
a  wedge.  The  distance  along  the  surface  (including  the  base),  S,  is 
normalized  by  the  top  surface  length,  L. 
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Figure  5.47:  Vibrational  temperature  jump  for  a  Mach  25  flow  of  nitrogen  about 
a  wedge.  The  distance  along  the  surface  (including  the  base),  S',  is 
normalized  by  the  top  surface  length,  L. 


201 


Table  5.6:  Computational  details  for  a  flow  of  nitrogen  about  a  wedge.  Total  CPU 
time  is  the  wall  time  multiplied  by  the  number  of  CPUs. 

Mach  10 


Krioo 

DSMC 

Cells 

Particles 

Time  Steps 

CPUs 

Total  CPU  Time* [hours 

0.002 

236,351 

57.5x10s 

360,000 

64 

6,483 

0.01 

215,090 

47.9  x10s 

250,000 

32 

1880 

0.05 

27,615 

9.5x10s 

200,000 

8 

248 

0.25 

7,137 

7.8  x10s 

200,000 

8 

200 

Krioc 

CFD 

Cells 

Iterations 

CPUs 

Total  CPU  Time*[hours 

0.002 

25,200 

30,000 

8 

130 

0.01 

25,200 

30,000 

8 

116 

0.05 

28,400 

30,000 

8 

126 

0.25 

30,175 

30,000 

8 

140 

Mach  25 


Klloo 

DSMC 

Cells 

Particles 

Time  Steps 

CPUs 

Total  CPU  Time' [hours 

0.002 

168,464 

31.7x10s 

350,000 

64 

3,696 

0.01 

120,985 

28.6x10s 

270,000 

32 

2,395 

0.05 

24,632 

4.9x10s 

200,000 

8 

144 

0.25 

6,041 

5.2x10s 

200,000 

4 

142 

Klloo 

CFD 

Cells 

Iterations 

CPUs 

Total  CPU  Time' [hours 

0.002 

29,050 

30,000 

8 

140 

0.01 

25,200 

30,000 

8 

138 

0.05 

28,400 

30,000 

8 

128 

0.25 

30,175 

30,000 

8 

132 

Approximate 
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Again,  there  are  no  significant  differences  between  the  nitrogen  and  argon  flows, 
despite  the  additional  presence  of  thermal  nonequilibrium  for  the  diatomic  gas. 


CHAPTER  VI 


Comparison  with  Experiment:  Hypersonic  Flow 

over  a  Flat  Plate 


6.1  Introduction 

The  previous  chapters  focused  purely  on  numerical  results,  with  CFD  simula¬ 
tions  being  compared  directly  to  DSMC  simulations.  In  particular,  the  walls  were 
assumed  to  be  fully  diffusive;  that  is,  the  gas  molecules  colliding  with  the  wall  were 
assumed  to  accommodate  fully  to  the  wall  conditions.  Hence,  an  accommodation 
coefficient  of  unity  was  used  for  the  CFD  slip  boundary  conditions.  In  this  chapter, 
two-dimensional  CFD  solutions  are  compared  with  experimental  measurements  of 
a  hypersonic  flow  of  nitrogen  over  a  flat  plate  [21],  Several  different  values  for  the 
accommodation  coefficient  are  evaluated.  In  addition,  the  CFD  solutions  are  also 
indirectly  compared  to  DSMC  solutions  of  the  same  flow  [66].  Thus,  the  relative 
accuracy  of  CFD  and  DSMC  can  be  evaluated  against  a  realistic  flow. 

6.2  Background  and  Experimental  Results 

The  experiment  was  conducted  at  the  University  of  Virginia  by  Cecil  and  Mc¬ 
Daniel  [21],  Measurements  of  the  flow  were  taken  using  planar  laser- induced  floures- 
cence  (PLIF).  A  hypersonic  flow  of  nitrogen,  with  a  Mach  number  of  approximately 
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Figure  6.1:  PLIF  image  of  hypersonic  flow  over  the  flat  plate  model  (from  Ref.  [21]). 

The  gas  expands  in  a  free  jet  from  the  left,  creating  a  barrel  shock  (near 
the  top  of  the  figure).  Note  the  light-colored  rays  extending  from  the 
orifice  upstream,  showing  the  radial  variation  of  the  velocity.  Note  also 
the  shock  attached  to  the  plate  leading  edge. 


11.9  at  the  leading  edge  of  the  flat  plate,  was  created  by  expanding  an  iodine-seeded 
flow  of  nitrogen  from  a  settling  chamber  through  a  thin,  circular  orifice  and  over  the 
model  in  a  continuously  evacuated  vacuum  chamber.  A  PLIF  image  of  the  flow  is 
seen  in  Figure  6.1.  Here,  the  gas  is  seen  expanding  in  a  free  jet  from  an  orifice  at 
the  left  of  the  image  to  the  right  over  the  model.  The  process  of  expanding  the  flow 
through  the  orifice  created  a  barrel  shock  (the  top  portion  of  which  is  clearly  seen 
at  the  top  of  the  figure)  terminated  by  a  normal  shock  (not  shown).  The  expanding 
flow  from  the  orifice  within  the  barrel  shock,  where  the  model  was  located,  was  not 
uniform,  but  instead  varied  radially.  On  the  left  of  the  figure,  the  light-colored  rays 
extending  from  the  orifice  upstream  of  the  model  shows  the  radial  variation  of  the 
velocity.  The  shock  formed  by  the  flow  over  the  model  is  also  clearly  shown. 

The  temperature  of  the  settling  chamber  was  approximately  300  K.  Expansion 
through  the  orifice  reduced  the  temperature  to  about  11.5  K  near  the  leading  edge 
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of  the  model.  This  temperature  is  sufficiently  low  that  no  vibrational  activation  is 
expected. 

The  density  at  the  leading  edge  was  about  3  x  10”4  kg/ m3,  giving  a  global  Knudsen 
number  (based  on  a  hard-sphere  mean  free  path  and  the  flat  plate  length  of  20mm) 
of  about  0.009.  This  is  comparable  to  the  global  Knudsen  number  value  of  0.01  for 
the  cases  discussed  in  previous  chapters.  In  those  cases,  CFD  gave  decent  results 
for  the  wedge  flow,  with  the  exception  of  the  region  near  the  leading  edge  where 
nonequilibrium  effects  were  significant. 

Velocity  magnitude  contours  and  streamlines  from  the  experimental  results  are 
shown  in  Figure  6.2.  Note  that  the  length  coordinates,  x  and  y,  are  nondimensional- 
ized  by  the  diameter,  cl  =  0.5  mm,  of  a  small  nozzle  built  into  the  plate  at  x/d  =  30, 
or,  equivalently,  x  —  15  mm.  In  further  experiments,  a  jet  was  issued  from  the 
nozzle  to  simulate  a  reaction  control  system.  The  flat  plate  leading  edge  is  located 
at  x/d  =  0,  and  the  trailing  edge  is  located  at  x/d  =  40.  The  streamlines  in  the 
freestream  again  illustrate  the  nonuniform  flow  due  to  the  source-like  nature  of  the 
flow  expanding  from  the  orifice. 

It  should  also  be  noted  that  there  is  a  pocket  of  very  low  velocity  gas  near  the 
trailing  edge  of  the  model  (between  x/d  =  40  and  50).  There  is  evidence  of  an  adverse 
pressure  gradient  at  this  location  that  will  be  discussed  in  further  detail  below. 

Experimental  data  for  the  time-averaged  velocity  was  provided  for  y/d  =  0  to 
15  at  several  locations  along  the  plate  from  x/d  =  0  to  40  in  increments  of  5.  This 
corresponds  to  y  =  0  to  7.5  mm  and  x  =  0  to  20  mm  in  increments  of  2.5  mm.  The 
uncertainties  of  the  U  and  V-velocity  components  were  estimated  to  be  on  the  order 
of  50  m/s  and  30  m/s,  respectively  [21]. 
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Conditions:  To=  300  K  Po=  1360Torr  Pback=  260  mTorr  d  =  0.5mm 


V  [M/SEC] 
800 
750 
700 
650 
600 
J  550 
■  500 
|  450 


Figure  6.2:  Measured  velocity  contours  and  streamlines  on  the  centerplane  of  the 
hypersonic  flow  over  a  flat  plate  (from  Ref.  [21]).  A  small  nozzle  is  built 
into  the  plate  at  x/d  =  30.  Note  the  radial  variation  in  velocity  upstream 
of  the  flat  plate,  and  the  pocket  of  low  velocity  gas  near  the  trailing  edge 
at  x/d  =  40. 


6.3  Computational  Results:  DSMC 

Padilla  and  Boyd [66]  conducted  DSMC  simulations  of  this  flow  using  pure  ni¬ 
trogen  and  the  variable  soft-sphere  (VSS)  DSMC  model  [4],  The  nonuniform  inflow 
velocity  profiles  used  for  the  DSMC  case  were  taken  from  the  experimental  data.  The 
temperature  and  the  particle  number  density  were  calculated  assuming  an  isentropic 
expansion  from  the  settling  chamber  conditions  of  T  =  300  K  and  p  =  1.79  atm. 
These  inflow  boundary  conditions  are  shown  in  Figure  6.3,  where  the  streamwise, 
or  ^-direction,  velocity  component  is  U  and  the  wall-normal,  or  ^/-direction,  velocity 
component  is  V.  Note  that  the  particle  number  density  has  been  converted  to  the 
mass  density.  The  variation  of  the  R-velocity  is  particularly  important;  while  the 
[/-velocity  varies  by  only  6  m/s  across  the  inflow  boundary,  the  R-velocity  varies  by 
nearly  150  m/s  (from  about  20  m/s  at  y  =  0  mm  to  about  160  m/s  at  y  =  7.5  mm. 
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V  [m/s]  Mass  Density  [1 0 4  kg/m3] 


(a)  Velocity  components  (b)  Temperature  and  mass  density 

Figure  6.3:  Inflow  boundary  conditions  for  a  hypersonic  flow  over  a  flat  plate.  The 
inflow  boundary  conditions  for  the  two-dimensional  simulations  vary  with 
vertical  distance  from  the  plate  (y)  to  match  the  radial  variation  of  the 
experimental  conditions. 


The  inflow  boundary  was  set  2  mm  upstream  of  the  flat  plate  leading  edge.  A  wall 
temperature  of  300  K  was  assumed. 

DSMC  simulations  were  computed  with  two  different  gas-surface  interaction  mod¬ 
els  and  tangential  momentum  accommodation  coefficients  of  0  to  1.0.  The  DSMC 
simulations  show  that  a  tangential  momentum  accommodation  coefficient  of  1.0  gives 
the  best  agreement  with  the  measured  data  for  the  streamwise-component  of  velocity 
([/),  while  a  tangential  momentum  accommodation  coefficient  of  0.75  gives  the  best 
overall  agreement  with  the  wall-normal  component  of  velocity  ( V ).  An  average  ac¬ 
commodation  coefficient  of  0.875  gave  good  overall  agreement  with  the  experimental 
data. 


6.4  Computational  Results:  CFD 


CFD  solutions  are  obtained  using  the  Gokgen  [29]  slip  boundary  conditions.  Sim¬ 
ulations  performed  using  the  Maxwell  and  Lockerby’s  wall-function  [49]  slip  bound- 
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ary  conditions  show  very  similar  results,  with  the  exception  that  the  Gokgen  slip 
conditions  tend  to  show  a  larger  slip  velocity.  This  is  not  unexpected;  previous  sim¬ 
ulations  for  similar  flow  conditions  (see,  for  example,  Figures  4.54(b)  and  4.55(b)) 
showed  that  the  velocities  predicted  by  all  three  slip  conditions  are  very  similar  away 
from  the  wall,  with  the  Gokgen  simulations  predicting  the  highest  slip  velocity.  The 
higher  slip  velocity  values  predicted  by  the  Gokgen  slip  conditions  agree  best  with 
DSMC  and  the  experimental  data,  for  these  flow  conditions,  as  will  be  shown  below. 

The  VHS  viscosity,  with  the  values  for  nitrogen  given  in  Table  3.1,  are  also 
used  here.  Although  the  DSMC  simulations  used  the  VSS  model,  there  are  no 
significant  differences  expected  at  the  low  temperatures  encountered.  In  addition,  a 
test  simulation  using  Blottner’s  curve  fit  for  the  viscosity  [7]  yielded  no  discernable 
difference  in  the  velocity  profiles. 

The  method  used  here  for  producing  mesh-independent  results  is  somewhat  differ¬ 
ent  than  that  used  previously.  Earlier  studies  were  mainly  concerned  with  obtaining 
correct  values  for  the  surface  properties.  Therefore,  the  surface  property  values  were 
used  to  determine  when  a  mesh-independent  solution  was  obtained.  It  was  found 
that  mesh-independent  solutions  were  most  sensitive  to  the  node  spacing  near  the 
wall  (for  both  the  wedge  and  the  cylinder  cases)  and  in  the  streamwise  direction  (for 
the  wedge  case).  The  objective  here,  however,  is  to  compare  the  flow  field  prop¬ 
erties  (specifically  the  velocity  components)  with  the  experimental  measurements. 
Therefore,  a  mesh-independent  solution  is  defined  as  one  for  which  successively  re¬ 
fined  meshes  produce  no  discernable  differences  in  the  velocity  profiles.  For  this  flat 
plate  flow,  the  flow  field  properties  (especially  the  V -component  of  velocity)  is  very 
sensitive  to  the  node  spacing  in  the  wall-normal  direction  in  the  entire  flow  field, 
rather  than  simply  near  the  wall.  There  is  also  some  sensitivity  to  the  node  spacing 
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in  the  streamwise  direction,  especially  near  the  leading  edge  where  the  gradients  are 
largest. 

In  a  similar  manner  to  the  DSMC  simulations,  here  the  accommodation  coefficient 
is  varied  from  0.5  to  1.0.  The  computational  domain  consists  of  a  rectangular  area 
from  y  =  0  to  20  mm  and  from  x  =  —2  to  30  mm. 

6.4.1  Flow  Field 

The  contours  for  the  value  of  the  continuum  breakdown  parameter,  Kiigll,  as 
computed  from  the  CFD  solutions  for  two  different  accommodation  coefficient  values 
(a  =  1.0  and  0.5),  are  shown  in  Figure  6.4,  and  the  density  contours  are  shown  in 
Figure  6.5.  It  is  apparent  that  the  flow  is  dominated  by  viscosity  effects;  in  other 
words,  the  flow  has  a  large  amount  of  nonequilibrium  present,  since  the  transfer  of 
momentum  due  to  velocity  gradients  is  an  inherent  nonequilibrium  process.  The 
largest  values  of  the  breakdown  parameter  are  found  in  the  shock  region,  and,  to 
a  lesser  extent,  the  wake.  It  is  interesting  to  note  that  the  value  of  the  breakdown 
parameter  exceeds  the  critical  value  of  0.05  in  most  of  the  flow  region  behind  the 
shock,  due  to  the  merging  of  the  shock  and  the  boundary  layer.  Nevertheless,  the 
CFD  solutions  do  give  good  agreement  with  the  experimental  data  in  these  areas,  as 
shown  below. 

Reducing  surface  accommodation  from  full  (a  =  1.0)  to  half  (a  =  0.5)  effectively 
reduces  the  effect  of  viscosity  on  the  flow  field.  For  full  accommodation  the  shock 
is  stronger,  and  the  compression  ratio  is  higher  (with  maximum  density  ratio  of 
about  2.0  compared  to  1.7  for  half  accommodation).  In  addition,  the  strong  viscous- 
shock  interaction  at  the  leading  edge  is  more  apparent  for  full  accommodation,  as 
can  be  seen  by  the  curvature  of  the  shock  in  that  area.  Cecil  and  McDaniel  [21] 
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X  [mm] 


(b)  a  =  0.5 

Figure  6.4:  Kiiqll  field  for  hypersonic  flow  over  a  flat  plate.  Significant  nonequilib¬ 
rium  effects  are  expected  at  the  leading  edge  and  in  the  wake.  Full  ac¬ 
commodation  increases  the  viscous  effects,  increasing  the  shock  strength 
and  the  shock-boundary  layer  interaction  that  leads  to  a  curved  shock  at 
the  leading  edge. 
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Figure  6.5:  Density  field  for  hypersonic  flow  over  a  flat  plate.  Full  accommodation 
increases  the  viscous  effects,  increasing  the  compression  across  the  shock, 
as  well  as  the  shock-boundary  layer  interaction  that  leads  to  a  curved 
shock  at  the  leading  edge. 
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mention  that  the  shock  wave  angle  of  18°  far  exceeds  the  value  of  4.8°  predicted 
by  inviscid  theory  and  that  an  effective  angular  displacement  of  13  —  14°  by  the 
boundary  layer  is  necessary  to  produce  the  18°  angle.  The  reduction  in  shock  angle 
as  the  accommodation  coefficient  is  reduced  also  demonstrates  the  reduced  effect  of 
viscosity  as  surface  accommodation  is  reduced. 

One  can  compare  the  velocity  magnitude  contour  plots  with  streamlines  from 
the  CFD  results  in  Figure  6.6  with  the  similar  plot  of  the  experimental  data  in 
Figure  6.2.  First,  note  that  the  nonuniform  inflow  conditions  appear  to  match  the 
experimental  results  quite  well.  Also  note  that  the  surface  accommodation  in  the 
experiment  appears  to  be  lower  than  1.0,  as  the  streamlines  near  the  leading  edge 
in  the  experimental  results  are  not  displaced  upward  as  much  as  is  shown  in  the 
numerical  results  with  full  accommodation.  Finally,  note  that  the  computational 
results  do  not  show  the  pocket  of  low  velocity  at  the  trailing  edge  that  is  apparent  in 
the  experimental  results.  (Recall  that  x/d  =  40  is  the  location  of  the  trailing  edge, 
and  is  equivalent  to  x  =  20  mm.)  It  appears  that  whatever  causes  the  low  velocity 
at  the  trailing  edge  (perhaps  a  region  of  high  pressure)  also  causes  the  streamlines 
to  be  displaced  upwards  much  more  in  the  experimental  results  than  is  shown  in  the 
numerical  results. 

6.4.2  Velocity  Comparisons 

The  flow  velocity  component  results  from  each  numerical  simulation  are  compared 
with  the  experimental  results  at  several  locations  along  the  flat  plate.  The  [/-velocity 
profiles  are  compared  in  Figure  6.7  and  the  V-velocity  profiles  are  compared  in  Figure 
6.8.  The  U  profiles  show  typical  boundary  layer  behavior,  including  the  non-zero 
velocity  slip  at  the  wall  due  to  the  rarefied  nature  of  the  flow.  The  V  profiles  include 
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Figure  6.6:  Computed  velocity  magnitude  contours  and  streamlines  for  hypersonic 
flow  over  a  flat  plate.  The  half-accommodation  results  compare  more 
favorably  with  the  experimental  results  in  Figure  6.2. 


214 


the  effects  due  to  the  shock. 

The  first  profiles,  at  x  =  0  mm,  show  good  agreement  between  the  experiment 
and  the  simulations.  Although  there  is  slightly  worse  agreement  in  the  V- velocity 
profile  as  y  increases,  the  numerical  results  remain  well  within  the  estimated  level  of 
uncertainty  in  the  experimental  data.  The  disagreement  as  y  increases  is  most  likely 
dne  to  incomplete  inflow  boundary  condition  modeling.  The  values  used  were  derived 
from  the  experimental  data,  which  is  limited  to  values  below  y  —  7.5  mm.  Although 
it  can  be  assumed,  by  looking  at  the  experimental  streamlines  in  Figure  6.2,  that 
the  radial  component  of  velocity  continues  to  increase  as  y  increases,  the  boundary 
conditions  implemented  do  not  assume  an  increasing  value  of  V  as  y  increases  beyond 
y  =  7.5  mm  (see  Figure  6.8(a)).  Nevertheless,  the  boundary  conditions  used  in  the 
simulations  are  assumed  to  be  adequately  correct. 

Very  near  the  leading  edge,  at  x  =  1.5  mm,  the  flow  solution  with  a  =  0.5  shows 
the  best  agreement  with  the  experiment,  for  both  U  and  V.  Further  down  the  plate, 
however,  up  to  x  —  12.5  mm,  flow  solutions  with  more  surface  accommodation, 
a  =  60  or  70,  agree  best  for  U,  while  the  solutions  with  a  =  50  continue  to  agree 
best  for  V.  Starting  at  x  =  15  mm,  the  agreement  for  U  progressively  worsens,  and 
the  experimental  U  profiles  at  x  =  17.5  and  20  mm  strongly  imply  the  presence  of  an 
adverse  pressure  gradient.  The  specific  cause  of  the  pressure  gradient  is  not  present 
in  the  simulations,  and  it  is  not  surprising  to  see  disagreement  in  the  flow  velocity 
predictions  near  the  trailing  edge. 

It  is  not  surprising  that  the  lower  surface  accommodation  values  agree  better 
nearer  to  the  leading  edge;  less  surface  accommodation  reduces  the  affect  of  viscosity 
(and  increases  the  velocity  slip) — it  was  shown  previously  for  the  wedge  that  CFD 
tends  to  underpredict  the  value  of  the  velocity  slip  near  the  leading  edge. 


Y  [mm]  Y  [mm] 
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■■■■  Experiment 
-  CFD  (o  =  0.50) 
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Figure  6.7:  Velocity  component  parallel  to  the  surface  (U)  for  a  hypersonic  flow 
over  a  flat  plate.  All  CFD  results  match  the  experimental  values  at 
the  leading  edge,  confirming  that  the  inflow  boundary  conditions  are 
correct.  Simulations  with  higher  values  of  surface  accommodation  match 
the  experimental  results  best  near  the  leading  edge. 
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.  Experiment 

- CFD  (a  =  0.50) 


(f)  x  =  15.0  mm 


.  Experiment 

- CFD  (a  =  0.50) 


-  CFD  (a  =  0.75) 


(g)  x  =  17.5  mm  (h)  x  =  20.0  mm 

Figure  6.7:  Velocity  component  parallel  to  the  surface  (U)  for  a  hypersonic  flow  over 
a  flat  plate  (cont.).  Downstream  of  the  leading  edge,  CFD  results  with 
a  =  0.75  agree  best  with  experiment.  The  experimental  velocity  profiles 
near  the  trailing  edge  strongly  suggest  the  presence  of  an  adverse  pressure 
gradient. 
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.  Experiment 

- CFD  (o  =  0.50) 


-  CFD  (a  =  0.75) 


.  Experiment 

- CFD  (a  =  0.50) 


(a)  x  =  0  mm 


(b)  x  =  2.5  mm 


.  Experiment 

- CFD  (a  =  0.50) 

- CFD  (a  =  0.75) 


.  Experiment 

- CFD  (a  =  0.50) 


-  CFD  (a  =  0.75) 


(c)  x  =  5.0  mm  (d)  x  =  7.5  mm 

Figure  6.8:  Velocity  component  normal  to  the  surface  (V)  for  a  hypersonic  flow  over  a 
flat  plate.  The  inflow  boundary  conditions  match  the  experimental  data 
where  present  (below  y  =  0.75  mm).  The  computed  velocity  profiles 
through  the  shock  match  the  experiment  best  for  a  =  0.5. 


Y  [mm]  Y  [mm] 


218 


.  Experiment 

- CFD  (a  =  0.50) 

- CFD  (a  =  0.75) 


.  Experiment 
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(e)  x  =  12.5  mm 


(f)  x  =  15.0  mm 


.  Experiment 

- CFD  (a  =  0.50) 


-  CFD  (a  =  0.75) 


.  Experiment 

- CFD  (a  =  0.50) 


(g)  x  =  17.5  mm  (h)  x  =  20.0  mm 

Figure  6.8:  Velocity  component  normal  to  the  surface  (V)  for  a  hypersonic  flow  over 
a  flat  plate  (cont.).  The  computed  velocity  profiles  through  the  shock 
match  best  with  experiment  for  lower  values  of  surface  accommodation. 
Simulation  results  do  not  agree  well  with  the  experimental  values  near 
the  trailing  edge  due  to  an  adverse  pressure  gradient  in  the  experimental 
results. 
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(a)  Pressure  (b)  Shear  stress 

Figure  6.9:  Surface  pressure  and  shear  stress  for  a  hypersonic  flow  over  a  flat  plate. 

Decreasing  the  surface  accommodation  decreases  the  viscous  effects,  de¬ 
creasing  the  pressure  and  shear  stress. 

The  previous  DSMC  results  [66]  also  showed  that  simulations  employing  lower 
values  of  the  accommodation  coefficient  tended  to  agree  better  for  the  H-velocity. 
However,  the  DSMC  results  agreed  best  for  higher  levels  of  surface  accommodation 
(cr dsmc  —  1-0  rather  than  c Tcfd  =  70  for  U,  and  (Jdsmc  —  0-75  rather  than  (Tcfd  = 
0.5  for  V). 

6.4.3  Surface  Properties 

Although  there  are  no  experimental  results  for  surface  pressure  and  shear  stress, 
it  is  instructive  to  compare  the  numerical  results  at  the  different  values  of  accommo¬ 
dation  coefficient.  Figure  6.9  compares  the  surface  pressure  and  shear  stress  along 
the  flat  plate.  The  profiles  are  very  similar  to  those  for  the  wedge;  a  large  peak  at 
the  leading  edge  is  gradually  reduced  along  the  surface.  As  expected,  lower  levels  of 
surface  accommodation  reduce  the  shear  stress  on  the  surface.  The  surface  pressure 
is  also  reduced  with  lower  values  of  surface  accommodation  as  the  pressure  on  a  flat 
plate  is  due  solely  to  viscous  effects. 
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Figure  6.10:  Velocity  slip  for  a  hypersonic  flow  over  a  flat  plate.  Reduction  in  sur¬ 
face  accommodation  increases  velocity  slip.  The  large  decrease  in  the 
experimental  velocity  slip  near  the  trailing  edge  is  most  likely  due  to  an 
adverse  pressure  gradient. 


The  experimental  values  are  extrapolated  to  the  surface  and  plotted  with  the 
computational  results  for  velocity  slip  in  Figure  6.10.  As  the  flow  progresses  along 
the  plate  the  large  value  at  the  leading  edge  is  reduced  to  an  almost  constant  value. 
The  large  reduction  in  experimental  velocity  slip  near  the  trailing  edge  is  also  due 
to  the  adverse  pressure  gradient.  Extrapolating  the  experimental  velocity  slip  from 
the  U  profiles  show  that  near  the  middle  of  the  plate,  at  about  x  —  10  mm,  the 
experimental  values  fall  between  the  flow  solutions  with  a  =  0.75  and  cr  =  0.5. 

6.5  Computational  Details 

The  computational  details  for  the  simulations  of  a  hypersonic  flow  of  argon  about 
a  wedge  discussed  in  this  section  are  shown  in  Table  6.1. 
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Table  6.1:  Computational  details  for  a  hypersonic  flow  over  a  flat  plate.  Total  CPU 
time  is  the  wall  time  multiplied  by  the  number  of  CPUs. 


Cells 

Iterations 

CPUs 

Total  CPU  Time* [hours] 

78,375 

6,000 

8 

36 

Approximate 


6.6  Summary — Hypersonic  Flow  over  a  Flat  Plate 

CFD  simulations  of  a  hypersonic  flow  over  a  flat  plate  are  compared  with  ex¬ 
perimental  results  and  previously  obtained  DSMC  results.  Although  there  is  some 
evidence  of  flow  separation  near  the  trailing  edge  of  the  plate  in  the  experimental 
results,  CFD  velocity  data  agrees  fairly  well  with  the  experimental  data  at  other 
locations  along  the  plate.  For  the  velocity  component  parallel  to  the  flow  (£/),  re¬ 
sults  obtained  using  a  surface  accommodation  coefficient  value  of  0.5  agrees  better 
with  the  experimental  results  near  the  leading  edge,  while  further  along  the  plate  an 
accommodation  coefficient  value  between  0.75  and  1.0  provides  best  agreement.  The 
lower  accommodation  value  near  the  leading  edge  effectively  increases  the  slip  veloc¬ 
ity;  CFD  also  under-predicted  the  velocity  slip  near  the  leading  edge  of  the  wedge. 
For  the  wall  normal  velocity  component  (D),  an  accommodation  coefficient  of  0.5 
gives  best  agreement  with  experimental  data.  Previous  DSMC  simulations  showed 
best  agreement  with  the  experimental  data  for  higher  values  of  the  accommodation 
coefficient. 

Less  surface  accommodation  significantly  decreases  the  predicted  pressure  and 
shear  stress  near  the  leading  edge. 


CHAPTER  VII 


Conclusions 


7.1  Summary 

The  primary  goal  of  this  study  was  to  quantify  the  difference  between  DSMC 
and  CFD  simulations  in  determining  the  effects  of  different  levels  of  nonequilibrium 
on  the  surface  properties  of  pressure,  shear  stress  and  heat  flux  of  a  body  under 
hypersonic  flow  conditions. 

This  detailed  computational  study  of  the  effects  of  nonequilibrium  on  the  surface 
properties  of  a  hypersonic  vehicle  considered  the  flow  about  a  cylinder  and  an  in¬ 
finitely  sharp  wedge.  Several  different  flow  regimes,  from  the  continuum  and  into  the 
transitional,  were  considered — free  stream  Knudsen  numbers  were  0.002,  0.01,  0.05 
and  0.25 — with  Mach  numbers  of  10  and  25.  The  effects  of  translational  nonequil¬ 
ibrium  were  isolated  by  considering  flows  of  argon;  thermal  nonequilibrium  effects 
were  included  when  considering  flows  of  nitrogen  over  the  same  bodies  at  the  same 
flow  conditions.  Validation  of  the  CFD  code,  as  well  as  the  effect  of  different  levels 
of  surface  accommodation,  was  shown  by  considering  a  nitrogen  flow  over  a  flat  plate 
and  comparing  the  simulation  results  with  experimental  data. 

Comparison  of  CFD  and  DSMC  results  for  the  flow  of  argon  about  a  cylinder 
showed  that  the  surface  properties  of  pressure,  shear  stress  and  heat  transfer  rates 
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were  very  similar  for  the  lower  Knudsen  number  flows  where  the  continuum  hypoth¬ 
esis  is  valid,  as  expected,  while  the  results  diverged  in  the  higher  Knudsen  number 
cases.  The  surface  pressure  was  least  affected  by  continuum  breakdown,  as  quantified 
by  the  gradient-length  local  Knudsen  number,  among  those  properties  investigated, 
and  seemed  to  be  affected  only  by  continuum  breakdown  as  the  shock  and  boundary 
layer  merged  at  the  highest  Knudsen  number  flow.  The  shear  stress  was  most  influ¬ 
enced  by  nonequilibrium  effects.  The  addition  of  slip  velocity  and  temperature  jump 
boundary  conditions  greatly  improved  the  agreement  between  CFD  and  DSMC  at 
higher  Knudsen  numbers.  Several  different  types  of  slip  boundary  conditions  were 
examined,  and  the  best  agreement  appears  to  be  obtained  when  using  the  general¬ 
ized  slip  conditions  proposed  by  Gokgen  [29].  With  these  boundary  conditions,  the 
differences  in  total  drag  and  peak  heat  flux  predicted  by  CFD  and  DSMC  ranged 
from  less  than  1%  at  Kn  =  0.002  to  around  5%  at  Kn  =  0.25. 

For  the  case  of  a  simple  gas,  the  higher  velocities  associated  with  a  Mach  25  flow 
did  not  seem  to  increase  the  differences  between  the  CFD  and  DSMC  predictions. 
Although  the  extent  of  the  region  where  the  continuum  breakdown  parameter  ex¬ 
ceeded  the  critical  value  of  0.05  was  larger  at  the  higher  Mach  number,  the  predicted 
surface  properties  with  the  slip  boundary  conditions  still  remained  well  under  5% 
for  all  but  the  Kn  =  0.25  case,  where  the  peak  heat  transfer  rates  and  total  drag 
predictions  were  within  6%. 

Although  there  was  a  significant  amount  of  nonequilibrium  between  the  different 
thermal  modes  (translational,  rotational  and  vibrational)  when  considering  a  flow  of 
nitrogen,  the  trends  were  largely  similar  to  those  noted  when  considering  a  simple  gas 
with  no  internal  degrees  of  freedom.  The  pressure  and  shear  stress  were  least  sensitive 
to  the  nonequilibrium  effects,  while  the  heat  transfer  rate  was  most  sensitive.  Total 
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drag  differences  between  CFD  (with  the  best  slip  boundary  conditions)  and  DSMC 
remained  under  3%,  while  peak  heat  flux  differences  were  less  than  8%. 

For  the  cylinder  cases,  it  was  also  shown  that  as  the  Knudsen  number  increased, 
the  percentage  of  total  drag  due  to  skin  friction  forces  (versus  pressure)  increased 
as  well;  friction  drag  accounted  for  less  than  5%  of  the  total  drag  at  Kn  =  0.002  to 
about  20%  at  Kn  =  0.25.  Differences  in  drag  due  to  skin  friction  also  tended  to  be 
larger  than  differences  in  predicted  drag  due  to  pressure;  thus  the  larger  errors  at 
the  higher  Knudsen  numbers  were  due  mostly  to  errors  in  skin  friction  prediction. 

Differences  in  flow  property  prediction  were  generally  concentrated  within  about 
10  mean  free  paths  of  the  wall  surface,  or,  in  other  words,  within  the  Knudsen  layer. 

The  sharp-leading  edge  geometry  of  the  wedge  led  to  additional  flow  phenomena 
not  seen  with  the  cylinder,  which  in  turn  affected  the  surface  property  predictions. 

The  differences  in  total  drag  predicted  by  CFD  and  DSMC  for  the  wedge  were 
greater  than  with  the  cylinder  geometry,  due  to  the  larger  effect  of  the  friction  forces 
on  the  drag.  Additionally,  there  were  greater  differences  in  the  pressure  profiles  along 
the  wedge  surface,  but  the  effect  on  total  drag  was  relatively  small  (due  to  the  small 
angle  of  the  surface  with  the  flow).  Total  drag  differences  ranged  from  less  than  2% 
at  Kn  =  0.002  to  more  than  20%  at  Mach  10  and  34%  at  Mach  25  at  Kn  =  0.25. 
Unlike  the  cylinder  cases,  friction  forces  contributed  most  to  the  total  drag  of  the 
wedge.  For  Kn  =  0.002,  friction  forces  accounted  for  50%  of  the  total  drag,  which 
increased  to  nearly  90%  for  Kn  =  0.25.  Again,  differences  in  the  drag  due  to  friction 
was  much  higher  than  differences  in  the  drag  due  to  pressure. 

The  amount  of  nonequilibrium  near  the  leading-edge  significantly  affected  the 
prediction  of  temperature  gradients,  and  thus  there  were  significant  differences  in 
the  heat  transfer  rate  predictions — CFD  fails  to  adequately  predict  the  large  heat 
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fluxes  near  the  leading  edge  of  an  infinitely-sharp  wedge.  Peak  heat  flux  differences 
were  near  70%  (and  as  much  as  100%)  for  the  regimes  considered. 

As  with  the  cylinder,  there  were  no  significant  differences  between  the  nitrogen 
and  argon  flows  about  the  wedge  (although  drag  differences  were  slightly  less  with 
nitrogen,  as  with  the  cylinder  cases),  despite  the  additional  presence  of  thermal 
nonequilibrium  for  the  diatomic  gas. 

CFD  simulations  of  a  hypersonic  flow  over  a  flat  plate  were  compared  with  exper¬ 
imental  results  and  previously  obtained  DSMC  results.  Several  different  values  for 
the  surface  accommodation  coefficient  were  considered,  ranging  from  0.5  to  1.0.  Al¬ 
though  there  was  some  evidence  of  flow  separation  near  the  trailing  edge  of  the  plate 
in  the  experimental  results,  CFD  velocity  data  agreed  fairly  well  with  the  experi¬ 
mental  data  at  other  locations  along  the  plate.  For  the  velocity  component  parallel 
to  the  flow  ( U ),  results  obtained  using  a  surface  accommodation  coefficient  value  of 
0.5  agreed  better  with  the  experimental  results  near  the  leading  edge,  while  further 
along  the  plate  an  accommodation  coefficient  value  between  0.75  and  1.0  provided 
the  best  agreement.  The  lower  accommodation  value  near  the  leading  edge  effec¬ 
tively  increased  the  slip  velocity;  CFD  also  under-predicted  the  velocity  slip  near  the 
leading  edge  of  the  wedge.  For  the  wall  normal  velocity  component  ( V ),  an  accom¬ 
modation  coefficient  of  0.5  gave  best  agreement  with  experimental  data.  Previous 
DSMC  simulations  showed  best  agreement  with  the  experimental  data  for  higher 
values  of  the  accommodation  coefficient. 

Less  surface  accommodation  also  significantly  decreased  the  predicted  pressure 
and  shear  stress  near  the  leading  edge. 
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7.2  Contributions 

Although  there  have  been  some  studies  comparing  CFD  and  DSMC  simulations 
published  in  the  literature,  the  current  work  makes  several  unique  and  valuable  con¬ 
tributions  to  the  field  of  computational  simulation  of  hypersonic  aerothermo dynam¬ 
ics,  some  of  which  have  been  published  in  References  [50,  51,  52],  These  contributions 
are  enumerated  below. 

1.  Started  with  the  fundamentals.  Many  previous  studies  have  compared  CFD 
and  DSMC  with  experimental  data.  Thus,  they  have  included  complex  ther¬ 
mochemical  nonequilibrium  models  in  the  simulation.  The  present  study,  as  a 
purely  numerical  study,  has  focused  primarily  on  the  fundamentals  of  nonequil¬ 
ibrium  behavior  and  has  gradually  increased  the  complexity,  starting  with  a 
monatomic  gas,  argon,  and  progressing  to  a  diatomic  gas,  nitrogen.  The  effects 
of  each  type  of  nonequilibrium  on  the  surface  properties  are  then  quantified  as 
the  complexity  increases. 

2.  Studied  many  flow  regimes,  about  blunt  and  sharp  bodies.  Many  of  the  pub¬ 
lished  studies  are  limited  to  a  few  flow  regimes  or  body  geometries.  The  current 
work  is  more  comprehensive,  considering  flow  regimes  from  the  continuum  and 
into  the  transitional  regime  to  quantify  the  effects  of  the  degree  of  rarefaction; 
considering  two  different  flow  velocities  to  quantify  the  effects  of  larger  Mach 
number;  and  considering  two  types  of  geometry,  a  cylinder  and  a  wedge,  to 
quantify  differences  due  to  blunt-body  phenomena  versus  sharp  leading-edge 
phenomena. 

3.  Evaluated  the  effectiveness  of  several  types  of  CFD  slip  boundary  conditions 
and  compared  the  CFD  slip  values  with  the  DSMC  slip  values.  This  research 
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has  evaluated  the  effectiveness  of  several  CFD  slip  boundary  conditions,  in¬ 
cluding  one  only  recently  proposed  [49],  in  predicting  the  surface  properties 
of  a  hypersonic  vehicle.  The  actual  slip  quantities  (velocity  slip,  transla¬ 
tional/rotational  and  vibrational  temperature  jump)  predicted  by  these  bound¬ 
ary  conditions  are  also  compared  with  those  extracted  from  the  DSMC  simu¬ 
lations  for  each  flow  condition,  which  is  unique  to  this  dissertation. 

4.  Laid  the  foundation  for  further  studies  essential  to  the  design  of  hybrid  meth¬ 
ods.  Hybrid  methods  face  two  basic  problems;  determining  the  boundaries  be¬ 
tween  the  CFD  and  DSMC  domains  and  passing  information  from  one  domain 
to  the  other.  This  research  contributes  to  both  of  these  areas.  The  chosen  value 
for  the  continuum  breakdown  parameter’s  effectiveness  in  predicting  differences 
is  shown  by  comparing  the  breakdown  value  with  the  other  flow  properties.  An 
effective  hybrid  design  also  requires  that  the  different  submodels  used  in  both 
computational  methods  be  equivalent  as  much  as  possible;  thus  information 
passed  between  both  domains  is  as  equivalent  as  possible. 

5.  Showed  conclusively  that  flow  property  differences  near  the  surface  are  concen¬ 
trated  in  the  Knndsen  layer.  Unique  to  this  dissertation  are  the  results  that 
the  differences  between  CFD  and  DSMC  near  the  wall  are  concentrated  mainly 
in  the  Knndsen  layer,  defined  here  as  the  region  of  flow  10  mean  free  paths  or 
less  from  the  wall  surface. 

7.3  Future  Research 


The  present  research  is  the  most  complete  and  systematic  study  to  date  to  deter¬ 
mine  the  effects  of  continuum  breakdown  on  the  surface  properties  of  a  hypersonic 
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vehicle.  However,  there  are  still  additional  areas  in  which  the  work  can  be  expanded. 

Geometry 

The  first  area  in  which  to  extend  this  study  is  the  geometry  of  the  hypersonic 
body.  The  current  work  has  looked  at  two-dimensional  blunt-  and  sharp-leading  edge 
bodies,  specifically  a  cylinder  and  a  wedge.  It  was  seen  that  these  two  types  of  geom¬ 
etry  yield  large  differences  in  the  level  of  agreement  between  CFD  and  DSMC  predic¬ 
tions  of  the  surface  properties,  particularly  for  the  peak  heat  flux.  Multi-dimensional 
effects  can  be  explored  by  considering  axisymmetric  bodies  such  as  spheres  and  cones; 
fully  three-dimensional  bodies  representative  of  actual  hypersonic  vehicles  might  also 
yield  important  information. 

O’Brien  and  Lewis  [63]  have  proposed  the  use  of  power-law  shaped  leading-edges 
for  the  blunting  of  leading-edges  while  minimizing  shock  stand-off.  The  power-law 
shape  can  be  aerodynamically  sharp  or  blunt,  depending  on  the  actual  power-law 
exponent.  It  was  seen  in  this  dissertation  that  CFD  agrees  better  with  DSMC  for 
blunt-bodies,  with  larger  differences  for  infinitely  sharp  bodes.  Further  studies  should 
extend  this  understanding  to  power-law  shapes. 

Chemistry 

The  research  contained  in  this  dissertation  was  restricted  to  thermal  nonequil¬ 
ibrium  only  (translational,  and  rotational/vibrational).  At  the  higher  temperatures 
normally  encountered  in  hypersonic  flows,  chemical  reactions,  such  as  dissociation 
and  recombination,  become  important.  Further  research  in  this  area  would  require 
an  understanding  of  both  forward  and  backward  reactions  rates  and  the  methods 
with  which  they  are  modeled  in  CFD  and  DSMC.  These  submodels  would  need  to 
be  as  equivalent  as  possible.  The  effects  of  chemical  nonequilibrium  on  the  surface 
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properties  of  the  hypersonic  vehicle  could  then  be  studied. 

Studies  that  include  chemical  nonequilibrium  would  also  necessarily  consider  any 
linking  between  vibrational  energy  activation  and  dissociation  rates,  as  well  as  surface 
catalysis. 

Gas  Mixtures 

While  the  research  presented  here  has  been  limited  to  the  consideration  of  only 
simple  gases  consisting  of  one  species,  real  hypersonic  flows  typically  involve  multi¬ 
ple  species.  Consideration  of  multiple  species,  such  as  oxygen  and  nitrogen,  would 
require  the  additional  consideration  of  inter-species  diffusion  coefficients  as  well  as 
the  definitions  of  mixture  properties  versus  species  properties  in  both  computational 
methods. 

Flow  Conditions 

The  current  work  has  looked  at  flow  conditions  ranging  from  the  continuum 
into  the  transition  regime,  at  both  Mach  10  and  Mach  25.  It  has  been  seen  that 
the  vibrational  modes  are  not  largely  activated  for  the  Mach  10  flows,  and  that 
vibrational  nonequilibrium  does  not  affect  the  surface  properties  much  even  at  the 
higher  velocities.  Additional  work  might  consider  higher  velocity  flows,  similar  to 
the  Stardust  re-entry  conditions  [12],  where  vibrational  nonequilibrium  is  expected 
to  be  more  important.  Work  in  this  area  might  also  separate  the  contributions  of 
translational  and  vibrational  energy  modes  to  the  heat  flux  in  the  simulations  and 
compare  the  CFD  and  DSMC  predictions  of  each  component. 
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Rotational  Energy  Equation  in  CFD 

The  CFD  code  used  here,  LeMANS,  currently  includes  a  two-temperature  model, 
with  one  temperature  being  a  combination  of  the  translational  and  rotational  tem¬ 
peratures.  Work  is  currently  being  performed  to  implement  a  separate  rotational 
energy  equation  into  the  code.  It  would  be  instructive  to  determine  if  the  additional 
modeling  of  rotational  nonequilibrium  improves  the  agreement  between  CFD  and 
DSMC  for  the  nitrogen  flows. 

Variation  of  Surface  Accommodation 

Finally,  with  the  exception  of  the  flat  plate  flow  in  Chapter  VI,  the  current 
research  has  been  limited  to  considering  only  full  accommodation  at  the  surfaces, 
for  both  momentum  and  energy.  For  full  accommodation,  CFD  and  DSMC  tended 
to  agree  well  at  the  lower  Knudsen  numbers.  However,  it  was  seen  for  the  flat 
plate  flow  that  less  surface  accommodation  was  required  in  the  CFD  simulations 
to  achieve  good  agreement  with  the  experimental  data  than  was  required  in  the 
DSMC  simulations.  Additional  studies  of  the  flows  about  cylinders  and  wedges  might 
include  variation  of  the  amount  of  surface  accommodation,  and  perhaps  independent 
variation  of  the  momentum  and  energy  accommodation,  to  determine  any  additional 
effects  on  surface  properties. 
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ABSTRACT 


NONEQUILIBRIUM  HYPERSONIC  AEROTHERMODYNAMICS  USING  THE 
DIRECT  SIMULATION  MONTE  CARLO  AND  NAVIER-STOKES  MODELS 


by 

Andrew  J.  Lofthouse 


Chair:  Iain  D.  Boyd 


This  dissertation  presents  a  detailed,  computational  study  quantifying  the  effects 
of  nonequilibrium  on  the  surface  properties  of  a  hypersonic  vehicle  by  comparing 
Navier- Stokes-based  Computational  Fluid  Dynamics  (CFD)  and  direct  simulation 
Monte  Carlo  (DSMC)  simulation  results  for  the  flow  about  a  cylinder  and  a  wedge. 
Physical  submodels  contained  in  both  computational  methods  are  ensured  to  be 
as  equivalent  as  possible.  Translational  nonequilibrium  effects  are  isolated  by  con¬ 
sidering  a  monatomic  gas,  argon.  Thermal  nonequilibrium  effects  are  included  by 
considering  a  diatomic  gas,  nitrogen.  Several  different  flow  regimes  are  considered, 
from  the  continuum  into  the  transitional  (freestream  Knudsen  numbers  are  0.002, 
0.01,  0.05  and  0.25),  with  Mach  numbers  of  10  and  25.  Effects  on  surface  proper¬ 
ties  (total  drag  and  peak  heat  transfer  rate)  are  quantified  at  each  flow  condition. 
Flow  held  properties  are  also  compared.  Continuum  breakdown  parameter  values 
are  compared  with  other  how  and  surface  properties. 


The  effectiveness  of  several  types  of  CFD  slip  boundary  conditions  is  evaluated, 
and  the  velocity  slip  and  temperature  jump  (including  vibrational  temperature  jump) 
values  are  compared  with  those  extracted  from  DSMC  simulation  results.  The  slip 
conditions  of  Gokgen  (AIAA  Paper  1989-0461)  most  accurately  predict  surface  prop¬ 
erties,  while  the  slip  conditions  of  Lockerby  et  al.  ( AIAA  J.  43(6)  (June  2005), 
1391-1393)  agree  best  with  DSMC  slip  values. 

For  flows  of  argon  and  nitrogen  about  a  cylinder,  CFD  total  drag  predictions 
remain  within  6%  of  DSMC  predictions,  and  heat  flux  agreement  is  8%  or  better. 
For  flows  about  a  wedge,  total  drag  differences  range  between  2%  and  34%,  mostly 
due  to  friction  force  differences.  Peak  heating  differences  are  between  70%  and  100%; 
DSMC  predicts  a  much  higher  temperature  near  the  leading  edge  than  CFD. 

Flow  property  differences  near  the  wall  surface  are  shown  to  be  concentrated 
primarily  in  the  Knudsen  layer.  Validation  of  the  CFD  code,  as  well  as  the  effect  of 
various  levels  of  surface  accommodation,  are  shown  by  considering  a  nitrogen  flow 
over  a  flat  plate  and  comparing  the  simulation  results  with  experimental  data. 


